Spondyloarthritis, acute anterior uveitis, and Crohn’s disease have both shared and distinct gut microbiota

Patient and public involvement

Neither patients nor the public were involved in the design, conduct, reporting, or dissemination plans of our research.

Patient inclusion criteria

The German Spondyloarthritis Inception Cohort (GESPIC) is an ongoing prospective cohort initiated to study the course and long-term outcomes of SpA across its whole spectrum of clinical presentation, including but not limited to gut microbiome composition. In this study we cross-sectionally analyzed only the baseline data. Serum, stool, and peripheral blood mononuclear cell (PBMC) samples, demographic information, and clinical characteristics were collected.

Since September 2015, patients have been recruited in three different arms depending on their main condition: 1) established radiographic axSpA: patients were required to fulfill the modified New York criteria and be eligible to start biological disease-modifying anti-rheumatic drug (bDMARD) therapy, by presenting high disease activity (BASDAI ≥ 4 and/or ASDAS ≥ 2.1) despite previous treatment with nonsteroidal anti-inflammatory drugs; 2) Crohn’s disease¹⁹: patients were classified according to the Montreal classification including location and behavior of CD, and had been recently diagnosed; 3) acute anterior uveitis²⁰: patients with non-infectious AAU diagnosed by an ophthalmologist. CD and AAU patients were enrolled regardless of musculoskeletal symptoms and an experienced rheumatologist was responsible for the final diagnosis of SpA / no SpA for patients included in these cohorts. All cohorts were approved by the ethical committee (Charité-Universitätsmedizin Berlin, Berlin, Germany). All patients enrolled were at least 18 years of age and gave their written informed consent. Patients included in this analysis were naïve to or had not received treatment with bDMARDs for at least three months before the enrollment in the study, nor had they received systemic antibiotics for at least one month prior to their baseline stool sample. There were no other restrictions concerning therapy. Control individuals were selected from the Optiref study²¹, which consisted of patients with chronic back pain who went through a standardized rheumatologic examination where the diagnosis of SpA was ruled out. Individuals with CD, AAU or psoriasis were excluded from the control group.

16S rRNA amplicon sequencing

After storage at -80°C, fecal samples were defrosted on ice. Aliquots of fecal material (1ml) resuspended in RNALater were centrifuged and washed once with water to remove excess fixative and salt. Then, DNA was isolated using the ZymoBIOMICS DNA Miniprep Kit (Zymo Research) according to the manufacturer’s instructions. Bead Beating was performed four times for 5 minutes each. Amplification of the V4 region (F515/R806) of the 16S rRNA gene was performed according to previously described protocols²² across two separate experiments. Briefly, 25 ng of DNA were used per PCR reaction (30 μl). The PCR amplification was performed using Q5 polymerase (NEB Biolabs). The PCR conditions consisted of initial denaturation for 30s at 98°C, followed by 25 cycles (10s at 98°C, 20s at 55°C, and 20s at 72°C). Each sample was amplified in triplicates and subsequently pooled. After normalization PCR amplicons were sequenced on an Illumina MiSeq platform (PE300).

Taxonomic profiling and preprocessing

Computational analysis was carried out in R (v4.0.3). Demultiplexed reads were analyzed using the DADA2 pipeline (v1.25.2) and taxonomic assignment was based on the SILVA rRNA database (v138.1) at 80% bootstrapped confidence. A total of 2118 unique bacterial ASVs were identified across 435 samples, annotated at the genus- or family-level. The DECIPHER package was used for multiple sequence alignment and the phangorn package was used to build a phylogenetic tree. Samples failing to meet inclusion criteria, lacking important clinical metadata, or not passing quality control (i.e. library sizes <5000 and >38K total raw reads) were excluded from further analysis (SFig. 1A), resulting in 339 samples. Libraries were rarefied to 5000 reads to account for variability in sequencing depth, resulting in 2066 ASVs present in at least one sample. The phyloseq::tax_glom() function was implemented in a custom manner that binned ASVs at the lowest taxonomic annotation available, either genus (70% of ASVs) or family (remaining 30%), and separately according to phyla (100% of ASVs). After total-sum scaling (TSS) normalization, three datasets: (i) ASV relative abundances, (ii) genus- or family-level relative abundances, and (iii) phylum-level relative abundances, which could be further processed (i.e. filtered) as needed.

Alpha and beta diversity analysis and enterotyping

The phyloseq::estimate_richness() function was applied to raw sample counts in order to calculate the Shannon entropy as a measure of alpha diversity (within-sample variation). To assess beta diversity (between-sample taxonomic variation), the vegan::vegdist() function was applied to rarefied ASV profiles (after TSS normalization) to calculate pairwise Bray-Curtis dissimilarity matrices and the stats::cmdscale() function was used to perform a principal coordinates ordination analysis (PCoA). The DirichletMultinomial R package was used to enterotype samples from rarefied counts binned at the genus- or family-level, according to the procedure from Holmes, Harris, and Quince²³.

Differential abundance analysis

ASVs not present in more than 5% of samples or displaying mean relative abundances less than 1e-05 across all samples were filtered out, resulting in n=442 ASVs and n=123 binned higher-level taxa (annotated at either genus- or family-level, see Supplementary Table 1) for univariate statistical testing and confounder analysis. To examine a potential shared disease signature, all n=277 disease samples (“GESPIC”) were pooled and compared to the controls in an integrated differential abundance and association testing procedure. Samples were also grouped according to phenotype and phenotypes with n>40 samples per group (i.e. SpA-only, CD-only, AAU-only, SpA+AAU) were subjected to this procedure to resolve phenotype-specific associations. The phyloseq and ggtree packages were used to visualize the phylogenetic information of ASVs.

Custom linear models were built for all n=565 taxa using the lm function in base R with the formula relative_abundance ∼ disease_status * sequencing_experiment to account for (un)known technical variation. Model parameters were extracted, tested for significance using the default Wald method, adjusted for multiple testing (Benjamini-Hochberg procedure), and bootstrapped to obtain 95% confidence intervals using the parameters::model_parameters() function. As a robustness check, all taxa were tested for significant differential abundance using a blocked Wilcoxon permutation test (coin::wilcox_test() function with the formula relative_abundance ∼ disease_status | sequencing_experiment), and raw P values were adjusted for multiple testing.

Confounder and mediation analysis

Keeping the same sets of taxa and phenotypic groups described above, non-parametric effect sizes (Cliff’s delta and the Spearman correlation) were calculated for each possible pairing of taxonomic and clinical features (see Supplementary Table 2 for a complete list), then tested for significance independent of a disease signal using blocked Wilcoxon or Spearman tests from the coin package (formula relative_abundance ∼ covariate | disease_status), and adjusted for multiple testing (Benjamini-Hochberg procedure).

The mdt_simple() function from the JSmediation package was used as described in Yzerbyt et al.²⁴ to test specific causal mediation hypotheses. This approach used linear regression to estimate the paths (a, b, c, and c’) along a possible causal pathway involving three variables. The joint significance of these paths was assessed to determine if mediation was present, and if so what the direct and indirect contributions within the pathway were. To lower the burden of multiple testing, only configurations with significantly disease-associated taxa (as in Fig. 2) and covariates emerging from the above testing were considered.

Further linear models were built with variables tracking current and prior medication intake (see Supplementary Table 3, detailed therapies in Table 1 were not used) to estimate the impact of specific drug therapies on taxon abundances. Instead of using the case-control groupings above, samples were first grouped according to disease phenotype (groups with <40 samples were not analyzed), such that disease status was not a factor, and the formula relative_abundance ∼ intake_status * sequencing_experiment was used to test specific disease-drug differential abundance hypotheses (e.g. conventional synthetic (cs)DMARD use, which was only relevant in CD patients). In the group of all patients examining previous antibiotic use, the formula relative_abundance ∼ disease_status * sequencing_experiment + intake_status was used; in the control subset examining NSAID use, simple models (relative_abundance ∼ intake_status) were used as no other factors were relevant.

Clinical presentation of SpA, AAU, and CD cohorts

Demographic and clinical characteristics of our 277 patients (72 CD, 103 AAU, and 102 axSpA patients) and 62 back pain controls are detailed in Table 1. More than half of the AAU cohort and one fifth of the CD cohort presented with predominantly axial SpA, with just five cases from the CD and three from the AAU cohorts diagnosed as exclusively peripheral SpA. A large majority of AAU and axSpA cohort patients carried HLA-B27 while less than 10% of individuals in the CD and control groups did (Table 1 and SFig. 1B). Over 90% of all patients were naïve to biologic disease-modifying anti-rheumatic drugs (bDMARDs, or biologics). In terms of disease activity, only patients in the SpA cohort showed clearly high systemic inflammatory activity, with mean CRP levels of 14.0mg/l and Ankylosing Spondylitis Disease Activity Score (ASDAS) of 3.5±0.8. Patients in the CD cohort had relatively inactive disease with a mean Harvey-Bradshaw Index (HBI) of 3.3±3.9, and patients in the AAU cohort had mean CRP levels of 4.1mg/l, although 45 (43.7%) had an active episode of anterior uveitis at the time of enrollment.

High-level variation and taxonomic diversity of gut microbiota

We performed 16S rRNA sequencing and taxonomically profiled a total of 339 stool samples (277 patients and 62 disease-negative back pain controls). At the phylum level our disease cohorts were dominated by Firmicutes, followed by Bacteroidota, Actinobacteriota and Proteobacteria (SFig. 2A, 2B). At the genus level, we observed substantially more variation between individuals in all cohorts, with Bacteroides, Prevotella, and Faecalibacterium comprising the dominant genera in each cohort (SFig. 2A, 2C). Notably, back pain controls and patients from the axSpA and AAU cohorts had more Prevotella-dominant individuals on average than did the CD cohort, which had the highest proportion of Bacteroides-dominant individuals (SFig. 2B, 2C). Patients with CD phenotypes had the lowest alpha diversities (Fig. 1A, 1C), but beta diversity did not appear to correlate with disease phenotype (Fig. 1B) or technical variation due to the sequencing experiment (Fig. 1D).

• Download figure
• Open in new tab

Figure 1: CD phenotypes display lowest alpha diversity and beta diversity is correlated with the ratio of Prevotella:Bacteroides.

The Shannon entropy was selected to capture the alpha diversity and pairwise Bray-Curtis dissimilarities were used for a principal coordinates ordination analysis (PCoA) to evaluate beta diversity (see Methods). A) Alpha diversity grouped by primary cohort. B) Beta diversity, samples colored according to phenotype. C) Alpha diversity grouped by phenotype. D) Beta diversity, samples colored according to sequencing experiment with the spatial means of each group depicted (squares). E) Alpha diversity grouped by phenotype and enterotype (see Methods and SFig. 2). F) Beta diversity, samples colored according to the Prevotella-to-Bacteroides ratio²⁵ which correlated with the first projected ordination axis (Pearson’s r = -0.8, P<0.001), accounting for up to 10% of between-sample variation and indicating an inverse relationship between these taxa. AAU, acute anterior uveitis; CD, Crohn’s disease; (ax)SpA, (axial) spondyloarthritis.

• Download figure
• Open in new tab

Figure 2: Immune-mediated disease phenotypes share a depletion of Lachnospiraceae taxa and an SpA-driven enrichment in Collinsella.

Disease groups consisted of patients without concomitant disease unless explicitly stated (e.g. SpA+AAU); GESPIC refers to the group of all n=277 patients. Linear models were built for ASVs (n=442) and their genus- or family-level taxonomic bins (n=123), containing an interaction term to account for technical variation (see Methods). A) Phylogenetic relationships, regression coefficients, and prevalence of ASVs significantly associated with one or more disease states. Bolded taxa had an adjusted q<0.05 and an absolute coefficient estimate >0.001 (corresponding to a mean differential abundance between disease and control groups of .1%). The x-axis broadly represents evolutionary distance based on the 16S V4 amplicons. B) Regression coefficients for the bolded taxa in A, plus any bins which were significantly differentially abundant, colored by genus or family with 95% bootstrapped confidence intervals shown. Bolded taxa were significant using a blocked Wilcoxon test at a q<0.05 cutoff. C) Venn diagram containing the bolded taxa from B, summarizing the most robust associations. ASV, amplicon sequence variant; GESPIC, German Spondyloarthritis Inception Cohort; CD, Crohn’s disease; SpA, spondyloarthritis; AAU, acute anterior uveitis.

Taxonomic associations with immune-mediated inflammatory disease states

To uncover taxa potentially mediating the concomitance of SpA, CD, and AAU, we performed a differential abundance analysis with disease-negative controls. We observed significantly lower mean relative abundances of several Firmicutes in the disease group as a whole, mainly from ASVs belonging to the Lachnospiraceae family (Fig. 2A). Of these, Fusicatenibacter [ASV12], Anaerostipes [ASV21], and Blautia [ASV30] were among the most prevalent. One ASV from the Christensenellaceae family (a health-associated²⁶ taxon with known heritability²⁷) displayed the strongest overall depletion in the GESPIC group relative to the control group; when all seven ASVs in our data belonging to this family were binned, the magnitude of the depletion increased (Fig. 2B). Similar behavior was observed for Fusicatenibacter, Marvinbryantia, and Erysipelotrichaceae ASVs, but not for Bacteroides, Blautia, or Lachnospiraceae ASVs, which were, in contrast, not significantly differentially abundant at higher taxonomic levels. GESPIC patients also had significantly higher abundances of Collinsella and Flavonifractor, both of which tracked closely with single ASVs and were ultimately found to be driven by the SpA and CD phenotypes, respectively (Fig. 2C). In addition to Collinsella, the SpA group was enriched in several Lachnospiraceae ASVs, Subdoligranulum [ASV505], and Holdemanella [ASV121] (Fig. 2A), few of which were significant when binned (Fig. 2B). Similarly, the CD group was strongly enriched in a Ruminococcus ASV and strongly depleted in a Blautia ASV (Fig. 2A), neither of which were significant at their respective genus levels (Fig. 2B). In contrast, Phascolarctobacterium and Lachnoclostridium reached significance when binned, but not at the ASV level (in the SpA and CD groups, respectively). There were no ASVs or higher-level taxa significantly associated with the AAU phenotype and only two weakly enriched Lachnospiraceae ASVs in the SpA+AAU phenotype; however, adjusted mean differential abundances for these groups mostly followed the pattern of the other case-control comparisons (Fig. 2A). Taken together, our findings broadly reflect the known taxonomic diversity of the Lachnospiraceae family^28,29 at the clinical level, and reveal shared and disease-specific taxonomic signatures ranging from ASV to family resolution.

Impact of drug therapies in GESPIC patient microbiota

To clarify the disease-associated microbiota signals, we examined a wide range of clinically-relevant factors. We found no significant associations between any taxa and age, sex, or BMI, indicating our disease signals were unlikely to be confounded or mediated by demographic factors. Medication intake has been found to explain more variation in microbiome composition than disease status alone^16,17 and is therefore an important factor to consider in association studies. In contrast to much of the literature^30–33, no individuals in our study were taking bDMARDs at the time of sampling and a vast majority were naive to any biologic therapies, allowing us a better approximation of a true baseline disease signal. However, 70% (including back pain controls) were taking at least one glucocorticoid, csDMARD, or NSAID, and 22% were taking two concurrently (Supplementary Table 3).

As >95% of our SpA patients were taking NSAIDs, we could not examine the disease signal independent of this effect; however, 60% of our control group was taking NSAID monotherapy, which was associated with 1-2% more Fusicatenibacter [ASV12] and Subdoligranulum [ASV16] on average (Fig. 3A). Interestingly, these two genera specifically were found to predict good response to csDMARDs in rheumatoid arthritis patients³⁴. 38% of our SpA-only patients were additionally receiving proton pump inhibitors, which was associated with an increase in Phascolarctobacterium (Fig. 3B), likely confounding that signal (Fig. 2B). No other disease-associated taxa appeared to be sensitive to the drugs we were able to test.

• Download figure
• Open in new tab

Figure 3: Fusicatenibacter is enriched in NSAID monotherapy, and Akkermansia is enriched in CD patients previously treated with csDMARDs.

A) Linear regression coefficients for disease-associated taxa (Fig. 2B), estimating the differential abundance between medicated and unmedicated patients, either while adjusting for disease status and/or technical variation (two leftmost columns, see Methods), or resulting from a comparison within a specific subset of patients (specified underneath). Bolded taxa were most robustly associated with CD, SpA, or the shared signal (Fig. 2C). We had information tracking previous antibiotic use (>1 month and <3 months before sampling) for all individuals in our study. B) Regression coefficients from the same analyses in A, for any taxa with an absolute coefficient estimate >0.01 (corresponding to a mean difference between medicated and unmedicated groups of 1%), colored by genus or family with 95% bootstrapped confidence intervals shown. C) Alpha diversity among all n=262 patients for which information on current intake for all three drug therapies was available. ATB, antibiotic.

Nearly half our CD-only patients were on a form of csDMARD therapy at the time of sampling, which was associated with significantly less Prevotellaceae (Fig. 3B). When comparing csDMARD-naive CD patients to those that were currently or had previously received treatment, we found that treated individuals had increased Phascolarctobacterium and Akkermansia (Fig. 3B). Although species-level identification from amplicon data is a contentious practice we avoided here, it is worth noting that these individuals also had about 1% more of a Ruminococcus taxon [ASV17] identified as R. gnavus^28,35, enriched in our CD patients and elsewhere (using metagenomic data) in treatment-resistant UC patients³⁶, as well as SpA patients with concomitant CD³³.

Mediation of inflammation and disease activity by HLA-B27 and microbiota

Other clinically-relevant covariates included the presence or absence of the HLA-B27 antigen on immune cells, systemic inflammation (as measured by C-reactive protein, CRP), and disease activity scores (Fig. 4A). Fusicatenibacter was the taxa most strongly (negatively) correlated with CRP (Fig. 4B), and we estimated about 19% of the increase in CRP observed in our GESPIC group to be mediated by [the relative lack of] this taxon (Fig. 4C). Several Subdoligranulum, Roseburia, Lachnospira, and Faecalibacterium ASVs as well as the total (binned) Faecalibacterium abundance were significantly higher in HLA-B27+ individuals, whereas Escherichia-Shigella and Blautia [ASV83] abundances were significantly lower in those individuals (Fig. 4D).

• Download figure
• Open in new tab

Figure 4: Decreased Fusicatenibacter partially mediates increased CRP in patients with SpA or CD, and increased Faecalibacterium partially mediates higher disease activity in HLA-B27+ SpA patients.

A) Summary of Table 1 highlighting clinical factors differentiating our phenotypes. B) Non-parametric effect sizes and significance tests showing serum CRP correlations independent of disease status for all n=585 taxa across all n=339 individuals in our study. C) Mediation analysis where c represents the total effect of disease presence on serum CRP, and the sum of both direct (c’) and indirect (axb) paths, with the effect and significance of each component (a, b, c, c’) illustrated (see Methods). When accounting for (decreased) Fusicatenibacter abundances, serum CRP estimates were (significantly) further increased, from 6.7 mg/L to 8.28 mg/L in the GESPIC group compared to controls. *q<0.1, **q<0.01, ***q<0.001. D) Similar to B, with HLA-B27. E) HLA-B27+ individuals (mostly SpA and AAU) had increased inflammation and disease activity. F) Mediation analysis where c represents the total effect of HLA-B27 expression status on BASDAI scores and c’ the indirect effect accounting for the enriched Faecalibacterium [ASV4] observed in HLA-B27+. BASDAI, Bath Ankylosing Spondylitis Activity Index.

Since our disease phenotypes were significantly different in their CRP levels and expression of HLA-B27 (Fig. 4A), we hypothesized that some of our shared disease-associated taxa or HLA- and inflammation-associated taxa might serve as mediators. Regardless of group, none of the taxa we considered (nor summary features such as taxonomic richness and alpha diversity) were found to mediate the increase in CRP observed in HLA-B27+ individuals in our study (Fig. 4E). However, in the group of SpA-only patients and back pain controls (who had higher BASDAI scores on average than patients from the CD or AAU cohorts but few HLA-B27+ individuals), we found that Faecalibacterium [ASV4] mediated about 20% of the (modestly) significantly increased disease activity in HLA-B27+ individuals (Fig. 4F).

Pilot study of serum metabolites in immune-mediated pathologies

Lachnospiraceae taxa are among the main short-chain fatty acid (SCFA) producers in the gut^28,29, and dysbiosis also appears to be tightly correlated with the local balance between Tregs and pro-inflammatory Th17 cells³⁷, which respond to microbially-produced tryptophan metabolites³⁸. We performed a small, targeted pilot study to quantify these metabolites in the serum (see Supplementary Methods), and observed higher circulating levels of SCFA in GESPIC patients relative to controls, although not statistically significant (Fig. 5A). This was surprising but consistent with work that correlated increased circulating SCFAs with gut permeability and dysbiosis in diabetics³⁹. No metabolites correlated with disease-associated taxa, likely due to low statistical power. The kynurenine-to-tryptophan ratio (KYN/TRP, proposed as a metric of immune activation⁴⁰) was positively correlated with CRP in disease groups (Fig. 5B), in line with a previous study which found anti-TNFα therapy in SpA to be effective at reducing both clinical markers⁴¹. Serum serotonin concentration, found elsewhere⁴² to stratify CD patients into disease activity categories better than both CRP and KYN/TRP, was also strongly negatively correlated with BASDAI scores in individuals with SpA (Fig. 5C).

• Download figure
• Open in new tab

Figure 5: Serum serotonin and the kynurenine-to-tryptophan ratio correlate with disease activity scores and systemic inflammation.

A) Concentrations of circulating short-chain fatty acids (SCFAs, in μM/L), which were not significantly associated with any disease but appeared higher in disease patients (if available, median of the control group is shown as a dotted line). Some metabolites were below the limit of detection (<LOD), especially propionate, and thus not used for statistical analysis (see Supplementary Methods). Marginal density histograms include only valid measurements. B) Spearman correlations between systemic inflammation (CRP, in mg/L) and the kynurenine-to-tryptophan ratio (KYN/TRP), which were positive in disease cohorts, but not the control group (rho = -0.06). CRP values >5 mg/L (dotted line) implicate active disease. Both axes are square-root transformed to aid visualization. C) Spearman correlations between BASDAI scores and serum serotonin concentrations (in ng/ml). Overall Spearman correlation coefficients and P-values shown in B and C, blocked for groups shown in the respective legends. KYN, kynurenine; TRP, tryptophan.