Differences in Gut Bacteria Across Stages of Alzheimer's Disease in a Turkish Group

The study cohort consisted of 47 AD and 27 MCI (all amnestic) patients and 51 subjects who were nondemented controls (n = 125). To minimize the dietary confounding effect on the microbiome, we included healthy cohabiting spouses of the patients sharing the same diet as the controls. The control group therefore largely (n = 27) comprised partners of the patients. Participants were recruited in two health centers located in different cities. The cohort groups were statistically not different in terms of sex, but age and education factors were significantly different (Table 1) and therefore statistically adjusted in analyses. Expectedly, the groups were also different by cognitive tests, including the Mini-Mental State Exam (MMSE) and the Clinical Dementia Rating (CDR). Most AD participants had very mild or mild dementia, with CDR scores ranging from 0.5 to 3 (median CDR scores of 1 for AD, 0.5 for MCI, and 0 for the control group). The median MMSE scores were significantly higher in the control (score of 27) and MCI (score of 26) groups than the AD group (score of 16). A subset of AD patients (n = 12) was clinically asked to undergo lumbar puncture to ascertain diagnosis using CSF biomarkers, including Aβ42/Aβ40 ratio, phosphorylated tau (p-tau), and the p-tau/Aβ42 ratio (see Table S1 in the supplemental material). We collected medication information from the patient’s registry.

The gut microbiota was profiled using the V3-V4 hypervariable region of the 16S rRNA gene; The Nephele automatic pipeline denoised the paired-end sequences and assigned amplicon sequence variants (ASVs) according to the DADA2 algorithm (35). The Nephele pipeline produced both unrarefied and rarefied ASV tables. The rarefied table included a total of 3,486 ASVs in the table (10,769 sequences/sample) for downstream analyses.

We used the resulting rarefied ASV table’s relative abundances to plot taxonomic diversity. The phylum-level taxonomic analysis showed a typical human gut microbiota profile in terms of overabundance of Firmicutes, Bacteroidetes, and Proteobacteria (Fig. 1A). Together with Verrucomicrobia and Actinobacteria, the five phyla comprised 99% of all reads, but the phylum Proteobacteria was overrepresented in AD patient samples. Notably, the genus-level relative abundance distributions across samples showed Prevotella_9 and Bacteroides were the most abundant of the top 30 genera across the samples (Fig. 1B). To perform differential abundance analysis between samples, we sought concordance analysis among multiple tools. ANCOM-BC or ALDEx2, when using covariates in their models, both agreed that only Ruminoccus_unclassified is significantly differentially abundant among the groups (data not shown). Nevertheless, when we employed the limma-voom R package (with an age- and sex-adjusted false-discovery rate [FDR] of <0.05) we found that Prevotella_9, Bacteroides, and members of the Ruminococcaceae family were among the top most significantly differentially abundant taxa (ASVs) between the cohort groups (see the differentially abundant ASVs in Table S2, parts A and B, in the supplemental material); whereas these taxa did not reach significance when collapsed at the genus level (Table S2, parts C and D), pointing to the possibility that multiple species of these three taxa exist and may be associated with the health and disease states at the same time. A comprehensive comparative statistical assessment of multivariate and compositional methods (36) demonstrated ALDEx2 or alike tools suffer from low power, while limma-voom and Songbird in their own class were the best performers.

Alpha diversity indices (Shannon and inverse Simpson) (see Fig. S1A to D in the supplemental material) did not show significant differences after multiple testing corrections (Kruskal-Wallis, FDR > 0.05), but only the richness index, Chao1, showed a significant difference between MCI and the control group (pairwise Wilcoxon rank sum test, P = 0.008074).

We employed both relative abundance-based and recently developed compositionally aware tools, namely, DEICODE (37) and Songbird (38), to compare the composition and structure of bacterial communities in samples using multiple beta diversity indices (Bray-Curtis, Jaccard, and Aitchison). The principal-coordinate analysis (PCoA) showed separation of the three groups by both Bray-Curtis and Jaccard indices (Fig. 2A and B). We used the adonis2 function in the qiime2 plugin (q2-diversity) to perform permutational multivariate analysis of variance (PERMANOVA) with 999 permutations, and the included interaction terms (see Table S3 in the supplemental material) and separation of the groups were highly significant (P = 0.0001). Age and sex also significantly contributed to the total variance (P < 0.001), but the interaction terms were not significant. Furthermore, tests of dispersion between groups (PERMDISP) indicated only the dispersion MCI group is significantly heterogenous (pairwise comparisons, P = 0.033 for AD versus MCI, P = 0.024 for control versus MCI, and P = 0.672 for AD versus C), which may be attributed to unbalanced design. We added further support for the separation of the three groups representing the AD continuum from other ordinations. The canonical analysis of principal-coordinates (CAP) analysis unambiguously showed the three groups are distinct (trace statistic = 0.86855 and P = 0.001 for 999 permutations) (Fig. 2C). The final support in beta diversity was provided by the DEICODE analysis (robust Aitchison principal-component analysis [PCA]) (PERMANOVA, P = 0.02) (Fig. 2D), which indicated that the three groups are distinct, and the community clusters are largely driven by a subset of ASVs with the taxonomic assignments Prevotella_9, Bacteroides, an unclassified genus within Ruminococcaceae family (Ruminococcaceae_unclassified), and Escherichia/Shigella. Moreover, the co-occurrence analysis using SparCC showed that Prevotella_9 and Bacteroides were negatively correlated (correlation = −0.4445, FDR = 0.09355). Finally, the genus-level PCoAs showed partially overlapping clusters of these two taxa, while the groups overall were also significantly separated (PERMANOVA, P < 0.0001) (see Fig. S2A to C in the supplemental material). Collectively, the differential abundance and beta diversity analyses showed gut microbiota composition is associated with the disease status along the AD continuum. We therefore paid particular attention to these two taxa in the rest of the downstream analyses.

Enrichment analysis by multinomial regression embedded in the Songbird tool with regard to covariates (coded as formula: Age+Sex+Edu+MMSE+CDR+Groups(levels=(“C”, “MCI”, “AD”)), indicated that the natural log ratios of Prevotella_9 to Bacteroides and Prevotella_9 to Escherichia/Shigella significantly separated the AD group from the control group (Welch’s t test, FDR-adjusted P = 0.04), but not from the MCI group (Fig. 3A to D). (Note that the formula follows the Patsy formatting, as described in more detail in Materials and Methods.) Importantly, the Songbird tool excluded 25 samples from this analysis due to zero-rich abundances that do not allow for center-log ratio calculations. We therefore tested the natural log ratio of the top 25%, allowing us to include all samples in the analysis to the bottom 25% (sets 1 and 2 in Table S4, part A, in the supplemental material) of the ranked ASVs associated with the AD relative to the control group; also, the same ratios for the MCI group relative to the control group (sets 3 and 4 in Table S4, part B) and the ASVs enriched in each group were visualized with Qurro (39). Both sets of ranked log ratios revealed significant differences between the log ratios of features differentiating groups (Welch’s t test, FDR-adjusted P = 0.0002).

Considering the preceding results, we postulated that the gut microbiota profile along the AD continuum does not represent a single state, but rather distinct yet overlapping community types. We addressed this hypothesis using four unique methods: (i) partitioning around medoid (PAM)-based clustering (40), (ii) fitting Dirichlet multinomial mixture (DMM) models to partition microbial community profiles into a finite number of clusters (41) using the Laplace approximation, (iii) fitting latent Dirichlet allocation (LDA) (42, 43) using the perplexity measure, and (iv) analyzing topological futures of data density (44) based on the Mapper algorithm to capture subtle and nonlinear patterns of high-dimensional data sets and population-level stratification.

The PAM-based clustering identified three (k = 3) distinct clusters based on Gap statistics (see Fig. S3A in the supplemental material). PCoA of the sample abundances in the three clusters indicated significant separation of the clusters (PERMANOVA, P = 0.001) (Fig. 4). We confirmed the optimum number of clusters using both Jensen-Shannon and Bray-Curtis distance metrics (data not shown). The relative abundance of the genus Prevotella_9 dominated cluster 1, while the genus Bacteroides showed the highest relative abundances in the other two clusters (Fig. 4B).

Next, we employed Dirichlet multinomial mixture probabilistic community modeling using the DirichletMultinomial R package (41) and fitting genus-level absolute abundances. Based on Laplace approximation, three clusters (clusters 1, 2, and 3) represented the best model fit (Fig. 3B), which was congruent with the PAM-based clustering. The PCoA of these clusters and PERMANOVA pairwise tests further supported the existence of three distinct clusters within the microbial community (PERMANOVA, P = 0.01) (Fig. 4C). The genus Bacteroides was the most abundant taxon in the first two clusters, and the third cluster was dominated by Prevotella_9 (Fig. 4D). Notably, cluster 2 included a significantly higher abundance of Bacteroides (26.3%) than cluster 1 (9.9%) and cluster 3 (4.7%). In addition to highly enriched Bacteroides in cluster 2, the decreasing trend of Faecalibacterium abundance and elevated abundance of inflammation associated with Escherchia/Shigella suggested that cluster 2 can be named “Bacteroides2 (Bact2) enterotype,” as recently described (45, 46). Reportedly, the abundance of Bacteroides in Bact2 enterotype can reach as high as 78% in patients with inflammatory bowel disease and is associated with systemic inflammation. These results suggest that cluster 2 includes patients with aggravated systemic inflammation.

We also performed SIMPER analysis based on Bray-Curtis distance to identify taxa contributing most to dissimilarities between clusters (data not shown). Bacteroides, Prevotella_9, Faecalibacterium, and taxa within Ruminococcaceae family ranked among the top 10 taxa contributing the most to differences between the three DMM clusters. To examine which factors were associated with the DMM clusters, we analyzed the distribution of clinical metadata and diversity metrics within the clusters. Alpha diversity indices (Chao1, Shannon, and inverse Simpson) were statistically different between all three clusters after Benjamini-Hochberg FDR adjustment. In addition, the clinical factors (CDR, MMSE, age, sex, and education) were also all significantly different among the clusters (Fig. 5A to H), except age and MMSE, for which only cluster 2 reached significance (Kruskal-Wallis test followed by Dunn’s post hoc test, FDR < 0.05). (Fisher’s exact test was used for the parameter sex, which is not shown in the figure.)

As a third method, we tested the LDA potential to stratify gut microbiota of the cohort participants. This unsupervised machine learning technique is increasingly finding acceptance in the field of microbiome studies (47–49) for its unique ability to reveal latent or hidden groups within the data cloud. Figure S4 in the supplemental material shows the LDA model’s perplexity parameter and log-likelihood values to find the optimal number of clusters. Both parameters continued to partition the community without reaching a clear optimum; we therefore presented the first 10 subgroups, also to allow for comparison with the results of topological analysis. The perplexity parameter not reaching a clear optimum is unexpectedly consistent with recent publications using LDA in microbial ecology (47–49). Bacterial probability distributions (ranked by a probability of ≥1% in descending order) across the subgroups are displayed in Fig. 6A. Interestingly, of the 10 subgroups, two subgroups were dominated by Bacteroides (topics 1 and 5), and a subgroup (topic 2) was dominated by Prevotella_9, with 97% probability. These subgroups therefore resemble subgroups detected by PAM and DMM in terms of prevalence of Bacteroides and Prevotella_9. Unlike DMM and PAM, however, LDA detected a distinct subgroup (topic 10) with the top-ranking genus Escherichia/Shigella, which also included putatively opportunistic bacteria, such as Enterococcus and Klebsiella. Subgroups 4, 6, and 9 were conspicuous with the genera known to produce butyrate and acetate or were mucinophilic. Even though we present the first 10 subgroups (topics) here, we also examined higher-order subgroups and observed that the 10 subgroups are further partitioned into additional subgroups, such as subgroups in which the top-ranking probability of Lactobacillus and Akkermansia emerges. Finally, we plotted the Quetelet index by subgroups to infer associations between subgroups and the cohort groups (Fig. 6B). The Quetelet index estimates the relative change of the occurrence frequency of a latent subgroup among all the samples compared to that among the samples of the cohort groups. The index showed subgroups 1, 8, 9, and 10 are positively associated with the AD group. The subgroup 9 is enriched by the members of Ruminococcaceae family. The top-ranking Ruminococcaceae_UCG_002 and Akkermansia are more abundant in the AD group than the control group according to limma-voom analysis. Akkermansia overabundance in AD gut microbiota is counterintuitive but was previously reported by others (25), and this genus is more abundant in the gut microbiota of Parkinson’s patients as well (50). Subgroup 10, where Escherichia/Shigella is the top-ranking genus, is strongly associated with the AD group but negatively associated with other groups. Conversely, subgroups 2, 4, and 7, which are enriched by short-chain fatty acid (SCFA) producers, are positively associated with the control and MCI groups but negatively associated with AD.

Another and final method we employed to stratify gut microbiota was topological data analysis (TDA), based on the Mapper algorithm (51) embedded in the recently developed tmap tool (44). The tmap tool was developed for network representation for stratification and association study of high-dimensional microbiome data. After constructing TDA microbiome network using the Mapper algorithm (ordination, covering, and DBSCAN clustering), the workflow in the second step includes computation of a modified version of the spatial analysis of functional enrichment (SAFE) scores to map both the metadata and microbiome features into the TDA network to generate their vectors of SAFE scores. Vectors of SAFE scores are then used to perform ranking and ordination, as well as coenrichment relations to delineate relationship between metadata and microbiome features. To construct the TDA network, we first applied dimension reduction (filtering) in PCA using Bray-Curtis distance, followed by the above algorithm, and also repeated the entire analysis using Jensen-Shannon distance to reveal the effect of the distance metric, if any. To understand how driver taxa relate to each other and to the clinical metadata, we performed principal-component analysis (PCA) of SAFE scores. Figure 7A shows the TDA network and PCA (Bray-Curtis distance) of the taxa metadata based on SAFE scores (see Data Set S1, tab 1, in the supplemental material), respectively. We obtained a similar TDA network profile using Jensen-Shannon distance (Fig. 7B) and SAFE scores (Data Set S1, tab 2). The size of each marker is scaled according to its SAFE score, and only the top 30 bacterial species are shown in PCA figures for clarity. A node in the network represents a group of samples sharing similar bacterial genus profiles. Two given nodes are linked when common samples are shared between the two nodes. The TDA analysis using both distance indices resulted in very similar stratification profiles with the top 10 SAFE-scoring genera, including Prevotella_9, Bacteroides, Rumunococaceae_unclassified, species of Lachnospiraceae, and GCA90006675. Unsurprisingly, a few taxon rankings differed between the two profiles, such as Caprococcus_2 and Mollicutes_RF39_unclassified.

Furthermore, Fig. 8A and B show taxa and host covariates based on Bray-Curtis and Jensen-Shannon distances, respectively. Regardless of the distance metric, all three groups were clearly separated. The drivers of microbiome stratification (Prevotella_9, Bacteroides, and Ruminococcus_unclassified) are placed near the control, AD, and MCI groups, respectively, in both PCA figures. Of the clinical metadata, MMSE, sex, and education were grouped with the control group and coenriched with Prevotella_9, but also with the Prevotella_2, Haemophilus, and Lachnospiraceae_NK4B4 group. Conversely, CDR, age, and the AD group were clustered together and coenriched with taxa such as Subdoligranulum, Odoribacter, Bilophila, and Alistipes. The MCI group was coenriched with Ruminocoaceae_unclassified, Mollicutes_RF39_unclassified, Ruminocoaceae_UCG_005, and Lachnospiraceae_unclassified. However, some taxa, such as Odoribacter, were placed near the control group in the Jensen-Shannon distance PCA (Fig 8B), suggesting coenrichment of certain taxa can be somewhat influenced by the preferred distance metric.

To identify signature taxa, we constructed a random forest (RF) model based on selected features of gut microbiota and psychometric test scores (MMSE and CDR) that are typically used as a proxy in clinical diagnosis. We first used Songbird to select 300 ASV features (the top 25%) that differentiate between the healthy (control) group and the disease groups (MCI and AD). We then plotted the ASVs with the first 20 highest mean decrease in Gini values (Fig. 9A) and included ASVs with mean decrease in Gini values above the breakpoint curve in the RF analysis. We identified the following 9 ASVs above the breakpoint: Faecalibacterium (ASV45), Sutterella (ASV607), Coprobacter (ASV531), Bacteroides (ASV81), Anaerostipes (ASV364), Ruminoccocaceae_unclassified (ASV203), Lactobacillus (ASV65), Clostridium_sensu_stricto_1 (ASV118), and Ruminococcus_1 (ASV59). Notably, ASVs beyond the breakpoint line are largely the bacterial species responsible for the stratification of gut microbiota in samples such as Faecalibacterium, Bacteroides, and Ruminococcus_unclassified. We next calculated the diagnostic accuracy of the RF model by plotting the receiver operating characteristics curves (ROCs) for the above 9 taxa, MMSE, and CDR separately and in combination for each cohort group (Fig. 9B). The ROC values for these selected nine taxa were moderately accurate (area under the receiver operating characteristic curve [AUROC], 63%) (Fig. 9B), but when we included MMSE and/or CDR, we found that the RF model robustly classifies all three groups (groupwise AUROC range of 0.74 to 1.0) (Fig. 9B).

We used multivariate association with linear models (MaAsLin2) to assess the association between individual taxa and clinical metadata, including patients’ drugs (q ≤ 0.25). This analysis showed that Roseburia, Lactobacillus, and Fusicatenibacter were negatively associated with AD (see Fig. S5 in the supplemental material). Of the medication categories, there are several taxa found to be positively associated with antidepression drugs and statins. Blautia, Caprococcus, Butyricoccus, Dorea, Lachnospiraceae family members, and some members of Ruminoclostridium and Ruminococaceae, known to be butyrate producers, are all positively associated with antidepression drugs. Unexpectedly, we found that several taxa were significantly associated with statin medication, and of these taxa, Streptococcus and unclassified members of Erysipelotrichaceae were highly significantly associated with statin medication. We also observed the following taxa positively associated with statin medication: unclassified members of Ruminococaceae and Lachnospiraceae, Phascolarctobacterium, Desulfovibrio, Caprobacter, Bifidobacterium, Butyricoccus, Blautia, and Barnesiella.

The Istanbul Medipol University and Erciyes University Ethical Review Boards approved this study (approval no. 186/16.4.2015 and 85/20.02.2015, respectively). All participants were informed of the objectives of this study and signed a written consent form prior to their participation. The diagnoses of dementia and MCI due to AD were based on the criteria of the National Institute on Aging-Alzheimer's Association Workgroups on diagnostic guidelines for Alzheimer's disease (66, 67). Exclusion criteria for this study included history of substance abuse, any significant neurologic disease, and major psychiatric disorders, including major depression. Also, individuals were excluded who used commercial probiotics or antibiotics during the study period or within 1 month prior to providing a stool sample or who had undergone major gastrointestinal (GI) tract surgery in past 5 years. Both health centers followed the same protocols in recruiting cohorts and used kits from the same manufacturers to minimize the variations in wet lab procedures.

Cerebrospinal fluid (CSF) samples were included in the analyses from a subset of AD patients if the patient was requested to donate a CSF sample as part of the clinically mandated diagnostic protocol. CSF samples were collected in the morning after overnight fasting using spinal needles (22 gauge) and syringes at the L3/4 or L4/5 interspace. CSF was then aliquoted into 0.5-mL nonadsorbing polypropylene tubes and stored at −80°C until assay. Biomarker molecules in CSF (Aβ₄₂, phosphorylated tau (p-tau), and the p-tau/Aβ42 ratio) were measured consistent with the Alzheimer’s Association flowchart for lumbar puncture and CSF sample processing, and the biomarker levels were determined as previously described (68). Single 96-well enzyme-linked immunosorbent assay (ELISA) kits (Innogenetics, Ghent, Belgium) were used in quantitation. The cutoff values discriminating AD from controls were determined as reported before (69).

Stool samples from all participants were collected in the neurology clinics of the university training hospitals. The participants were given a collection kit included a sterile tube and provided a brief instruction for collection. Self-collected samples were placed within approximately 30 min of collection in −80 freezers and kept frozen until DNA extraction.

16S rRNA gene sequencing and PCR were performed as previously described (70), with minor modifications. Briefly, genomic DNA was extracted from 220-mg fecal samples using the QiaAmp DNA stool minikit (Qiagen, Germany) per the manufacturer's instructions, with the addition of bead beating (0.1-mm zirconium beads) and lysozyme and RNase A incubation steps.

To amplify the variable V3-V4 regions of the 16S rRNA gene, the primers 341 F (5′--3′) and 805 R (5′--3′) were used. MiSeq sequencing adaptor sequences were added to the 5′ ends of forward and reverse primers. Approximately 12.5 ng of purified DNA from each sample was used as a template for PCR amplification in a 25-μL reaction mixture by using 2× KAPA HiFi Hot Start ready mix (Kapa Biosystems, MA, USA). For PCR amplification, the following conditions were followed: denaturation at 95°C for 3 min, followed by 25 cycles of denaturation at 95°C for 30 s, annealing at 55°C for 30 s, and extension at 72°C for 30 s, with a final extension at 72°C for 5 min. Amplified PCR products were purified with the Agencourt AMPure XP purification system (Beckman Coulter), and Nextera PCR was performed by using sample-specific barcodes. The constructed Nextera libraries were then sequenced by the Illumina MiSeq platform using MiSeq reagent kit v2 chemistry. CCTACGGGNGGCWGCAG GACTACHVGGGTATCTAATCC

The paired-end 16S rRNA reads were first used cutadapt v1.9 program (71) for the process of quality filtering and trimming and uploaded on the DADA2 pipeline (35) integrated into the Nephele platform (72) (v2.0; http://nephele.niaid.nih.gov↗). Chimeric sequences are automatically removed by this pipeline, which generates both rarefied and unrarefied ASV abundance tables. We used the rarefied (10,769 reads/sample) ASV table in most downstream analysis due to large differences between some total sample reads, except for the scale-invariant DEICODE and Songbird. We removed any sequences that were classified as either being originated from Eukarya, Archaea, mitochondria, chloroplasts, or unknown kingdoms.

We included no-sample DNA extractions and no-template negative-control samples in every sequencing library prepared. Using reads in the negative-control samples as a reference, we identified and removed probable contaminant reads of 13 ASVs from the ASV table, as predicted by the Decontam R package (73) using the “prevalence” method. In this method, the binary coded features across samples are compared to the prevalence in negative controls to identify contaminants. Also, we sequenced the same amplicon of an AD sample three times to check the sequencing variation. Although both centers used same protocols and kits from the same manufacturer in sequencing, we sequenced amplicons amplified from two same genomic DNA templates again from AD samples at both centers to check the center-to-center sequencing concordance. No differences could be identified between the taxonomic compositions of the samples sequenced at both centers nor between the technical replicates (PCoA, PERMANOVA, P = 0.1).

Most numerical downstream analyses of ASV abundances were performed in the R environment (74). All P values, where appropriate, were adjusted for multiple testing using the Benjamini-Hochberg (false-discovery rate [FDR]) method. We measured within-sample microbial diversity (alpha diversity) using observed richness, Chao1, Shannon, and inverse Simpson in phyloseq (75) and tested using Kruskal-Wallis. To identify differentially abundant bacterial species, we employed the animalculus (59) and limma (76) R packages. We assessed microbial diversity between samples (beta diversity) by using multiple distance metrics, including Bray-Curtis, Jaccard, and canonical analysis of principal components (CAP). CAP analysis and the similarity percentages breakdown (SIMPER) procedure were performed using PRIMER v7 (77). Additionally, due to the compositional nature of the data, we also included robust Aitchison PCA, using the qiime2 DEICODE plugin (37) to calculate beta diversity with feature loadings. The resulting ordination was visualized using Emperor (78). We tested significance of beta diversity among groups using, again, qiime diversity plugin PERMANOVA.

Next, we used Songbird (38) for multinomial regression to rank species association with disease status with the following parameters: (formula: “MMSE+CDR+Sex+Edu+C(Group, Diff, levels=('C','MCI','AD'), –p-epochs 10000 –p-differential-prior 0.5 –p-summary-interval 1 –p-random-seed 3 –min-sample-count 1000 –min-feature-count 0). Of note, the formula structure follows Patsy formatting (https://patsy.readthedocs.io/en/latest/↗), such that Groups (C, MCI, AD) represent levels=[“healthy”, “mild”, “severe”] states, respectively. A null model was generated using the same parameters. The fitted model demonstrated better fit compared to the null model (pseudo Q² = 0.874027). Taxa ranks were visualized using Qurro (39). Significance was determined using a Welch’s t test between groups, performed by GraphPad Prism.

To identify microbial species associated with the clinical metadata, including patients’ medication, we performed multivariate association with linear models (MaAsLin2) (79). The control group was excluded from this analysis, as the members of this group were not normally prescribed these medications. We employed the R package MaAsLin 2.1.0 to perform per-feature tests. We log transformed relative abundances of microbial species and standardized continuous variables into Z-scores and binary encoded medication information before including them in the MaAsLin models (q < 0.25 for significance).

We employed clustering, probabilistic partitioning, and topological data analysis approaches for the stratification of gut microbiota in the samples. The partitioning around the medoid (PAM) approach (40) clusters samples by iteratively updating each cluster’s medoid. We assigned samples to community types using the function pam() in the R package “cluster” based on Bray Curtis and Jensen Shannon distances. The number of clusters was determined by Gap statistic evaluation. Departing from the clustering approach, we next used two distinct probabilistic methods to partition the microbiota landscape, namely, Dirichlet multinomial mixture (DMM) models (41) and latent Dirichlet allocation (LDA) (42, 43). Genus-level abundances were fitted to DMM models to partition microbial community profiles into a finite number of clusters, using the Laplace approximation as previously described (41, 80).

As a second probabilistic partitioning, we performed LDA, which is a multilevel hierarchical Bayesian model (42) otherwise used for collections of discrete data, such as text corpus analysis in linguistics. LDA is a generalization of Dirichlet multinomial mixture modeling in which biological samples are allowed to have fractional membership and distinct microbial communities have different microbial signatures. Thus, for each taxon, there is a vector of probabilities across all clusters that it can be assigned to. Each cluster, therefore, has a different probability of containing taxa, indicating the chance of microbes in a particular subgroup (strata) co-occurring due to community assembly dynamics. To fit the model, we used Gibb’s sampling with the R package MetaTopics (v1.0) (81). The relative abundances of a genus collapsed table with abundances of more than 0.1% and 5% sample prevalence were input to the model. We plotted perplexity measure and log-likelihood values to estimate model performance and the optimal number of topics (subgroups of microbial assemblages) using 5-fold cross-validation. However, we observed that both parameters continued to improve with increasing subgroup number without a clear optimum, except the first jump in perplexity was near 10 topics. We therefore picked the first 10 topics for the sake of interpretability.

The final method we applied was topological data analysis (TDA) based on the Mapper algorithm (51) and network representation for stratification and association of study of high dimensional microbiome data, all integrated into the tmap tool (44). The framework enables us to reveal the association of taxa or metadata within the entire network and to identify enrichment subnetworks of different association patterns. Conceptually, the Mapper algorithm transforms a distance matrix and represents the shape of the data cloud in an undirected network. A node in the network represents a group of samples with similar microbiome profiles, and if common samples between nodes are shared, then the nodes are linked. Next, a modified special analysis of functional enrichment (SAFE) algorithm maps the metadata and taxa into the network. Finally, vectors of SAFE scores can be used in ordination to rank the driver taxa and their relationship with the metadata.

To identify a microbial signature of severity of cognitive impairment in the AD continuum, we implemented a machine learning procedure. We first took advantage of the songbird tool to select features, including the covariates and healthy (control) and disease states (AD+MCI) in the model formula. We subsequently fit the list of ASVs selected this way into random forest models. We plotted the area under the receiver operating characteristic curve (AUROC) to show prediction performance of the models. To create the classifiers, a random forest constituted of 500 trees was computed using the default settings of the “randomForest” function implemented in the randomForest R package (v4.6-7). Mean decrease in Gini values were averaged for each ASV among the 100 random forest replicates. The ASVs with the first 20 highest mean decrease in Gini values were plotted. ASVs with mean decrease in Gini values above the breakpoint curve were chosen to be part of the classifier. Breakpoints were estimated using the “breakpoints” function included in the strucchange R package (52). We subsequently fit the list of ASVs selected this way with or without psychometric test values (i.e., MMSE and CDR) into random forest models and bootstrapped for 100 times. We plotted the area under the receiver operating characteristic curve (AUROC) to show prediction performance of the models.

The 16S rRNA sequences generated by this study have been submitted to the NCBI BioProject database (https://www.ncbi.nlm.nih.gov/bioproject/↗) under accession no. PRJNA734525↗.