Stools were collected from participants to an existing larger cohort study on brain aging in 18 memory clinics in Eastern Lombardy, Italy. The parent study aimed at assessing the added value of amyloid imaging in the clinical work-up of patients with cognitive complaints (the Incremental Diagnostic Value of Florbetapir Amyloid Imaging [INDIA-FBP] study) (Altomare et al., 2018). After completion of the INDIA-FBP procedures, 150 patients and controls were asked to further contribute with their samples of stools and blood. As previously reported (Cattaneo et al., 2017), stool samples were collected from subjects at their own home in a sterile plastic cup, stored at −20°C, and delivered to IRCCS Fatebenefratelli Institute in Brescia within the following 24 h, where they have been stored at −20°C until their processing. The stool samples used in the current study were from a subgroup of 40 subjects with cognitive performance from normal to dementia.

The study was approved by the Ethics Committee of “Comitato Etico dell’IRCCS San Giovanni di Dio – Fatebenefratelli” (Brescia, Italy) under registration number 57/2014. Written informed consent was obtained from all participants.

DNA was extracted from 180 to 200 mg of frozen stool using the QIAamp DNA Stool Mini Kit (Qiagen Retsch GmbH, Hanover, Germany) and according to the manufacturer’s instructions. Bead-beating homogenization by TissueLyser II (Qiagen Retsch GmbH, Hanover, Germany) was performed to mechanically disrupt fecal samples before DNA extraction. The samples were homogenized for 10 min at 30 Hz. DNA was then quantified using a NanoDrop ND-1000 spectrophotometer, and then stored at + 4°C for subsequent analyses. All the stool samples were processed at the same time.

Bacterial DNA was amplified and purified according to 16S Metagenomic Sequencing Library Preparation protocol by Illumina. The first step consists in the amplification of the regions V3 and V4 of the bacterial ribosomal RNA 16S gene, by using the suggested primers (Forward Primer = 5′TCGTCGGCAGCGTCAGATGTGTATAAGAGACA GCCTACGGGNGGCWGCAG; Reverse Primer = 5′GTCTCGT GGGCTCGGAGATGTGTATAAGAGAC AGGACTACHVGGG TATCTAATCC) and the suggested cycling conditions (3′ at 95°C; 25 cycles: 30′′ at 95°C, 30′′ at 55°C, 30′′ at 72°C; 5′ at 72°C). The amplicon DNA was immediately purified with a magnetic bead step, washed in 80% ethanol and resuspended in Tris-HCl 10 mM. Amplicon DNA integrity was assessed by 2.0% agarose gel electrophoresis on gels containing 0.8 mg/mL ethidium bromide and stored at −20°C. Within 1 week from storage, amplicon DNA was uniquely dual-indexed, by using the suggested indices (Nextera XT) and the suggested cycling conditions (3′ at 95°C; 8 cycles: 30′′ at 95°C, 30′′ at 55°C, 30′′ at 72°C; 5′ at 72°C). The resulting indexed DNA was immediately purified with a magnetic bead step, washed in 80% ethanol, resuspended in Tris-HCl 10 mM and stored at −20°C. Within 1 week from storage, indexed DNA was quantified using fluorometric quantification (Qubit, Invitrogen) and the amplicon length was determined by using a Bioanalyzer DNA 1000 chip (Agilent). DNA was then normalized to 4 nM, pooled, denatured with NaOH 0.1N, diluted to the final concentration of 10 pM and loaded into the MiSeq v3 cartridge (Illumina). A paired-end read of 300 cycles per read was performed. Samples from cognitively intact persons and demented patients were sequenced in two different runs. Data are publicly available at www.ebi.ac.uk/ena/data/view/PRJEB35434.

The paired-end MiSeq Illumina reads (2 × 300 bp) were processed by using four popular tools: QIIME2 (Bolyen et al., 2019) (version 2018.8), Bioconductor (Callahan et al., 2016b) (version 29 October 2018), USEARCH (Edgar, 2010) (version 11.0.667) and mothur (Schloss et al., 2009) (version 1.43.0). For QIIME2, Bioconductor and mothur, the 64 bit version was used while for USEARCH (Edgar, 2010), only a memory-confined 32-bit version is freely available for academic use. QIIME 2 was installed in a conda environment. All the pipelines were run on both a Linux workstation (Ubuntu 14.04.5 LTS) equipped with Intel CPU 8 × 3.70 GHz processors and 31.3 GB of RAM and a MacBook Pro (15-inch, 2018) with 2.6 GHz 6-core Intel Core i7 processor, 16 GB of RAM and macOS Mojave, version 10.14.6. In QIIME2 we followed the “Moving Pictures” tutorial (accessed 3 January 2019), in Bioconductor the pipeline reported in https://f1000research.com/articles/5-1492/v2 (version: 29 October 2018), in USEARCH (Edgar, 2010) the instruction of the online tutorial¹ (accessed 4 January 2019), and in mothur the SOPs in https://www.mothur.org/wiki/MiSeq_SOP (accessed 09/27/2018). Details for each pipeline are shown in Figure 1 and the list of commands used for each tool is reported in Supplementary Methods 1–6. Briefly, sequencing Illumina MiSeq data were already demultiplexed. The parameters used for the quality check of the 16S rRNA data were those suggested in the corresponding tutorials. In general, forward and reverse primers, reads containing ambiguous bases or homopolymers greater than eight base pairs in length as well as chimeras were removed in all pipelines. Moreover, in QIIME2 we set a maximum number of expected errors equal to 2 and reads truncation if the quality score was less than 2. In Bioconductor we used a maximum of 2 and 5 expected errors per-read for forward and reverse reads, respectively, and reads truncation if the quality score was less than 2; moreover, features with ambiguous phylum annotation were discarded. For both QIIME2 and Bioconductor, the software packages Divisive Amplicon Denoising Algorithm 2 (DADA2) (Callahan et al., 2016a) was used to infer true biological sequences from reads. Since the samples were sequenced in two different runs and that different runs may have different error profiles, the denoising process was applied separately to the two runs. In QIIME2, the “dada2” plugin uses the DADA2 R library and has been implemented to perform sequence denoise, dereplication, and chimeras filtering. All these steps correspond to individual commands in Bioconductor. For QIIME2, the pipeline included also the extraction of the representative sequences using “feature-table” and their classification by taxon using the “feature-classifier” (Supplementary Methods 1). For Bioconductor, the pipeline included also the taxonomic classification of the ASVs using “assignTaxonomy,” the construction of the phylogenetic tree using the phangorn R package (v2.4.0) (Schliep, 2011) and the creation of a single data object that can be used for graphical and statistical purposes using phyloseq R package (v1.24.2) (McMurdie and Holmes, 2013) (Supplementary Methods 2). The UPARSE (Edgar, 2013) and the UNOISE (Edgar, 2016) pipelines are both implemented in USEARCH and have been used to defined OTUs with a similarity threshold of 97 and 99%, respectively, as recommended by the author of USEARCH². All the steps preceding the clustering, namely merging, filtering (where a maximum number of expected errors equal to 1 was set) and dereplication, are the same between the two pipelines (Supplementary Methods 3). In mothur, sequences were merged using “make.contig” command, quality filtered by “screen.seqs” command, dereplicated, aligned to the reference database SILVA v132 (Quast et al., 2013) and pre-clustered allowing 1 difference for every 100 bp of sequence (four differences in our case) (Supplementary Methods 4). After removing the undesirables features (i.e., chloroplasts, mitochondria, archaea, eukaryota, or unknown) using the “remove.lineage” command and chimera by applying the “chimera.vsearch” command, unique sequences were clustered into OTUs at 97 and 99% identity.

SILVA (version 132) (Quast et al., 2013)³ and RDP (version 16) reference databases were customized following the instructions on the respective tutorials for QIIME2 and mothur (Supplementary Methods 5, 6) and were downloaded for Bioconductor (SILVA v132⁴, RDP v16⁵). As SILVA version 132 was not available for UPARSE, taxonomic assignment of OTUs was performed on mothur (Supplementary Method 3), where bacterial classification identity was set to 80%. ASVs and OTUs found in a single sample (also named as singletons) were discarded.

In order to facilitate the data analysis and reproducibility, the list of all commands used for each pipeline and the number of reads assigned to taxa, are included as Supplementary Methods 1–6 and Supplementary Table 1.

All statistical analyses and graphs were performed using GraphPad Prism Software (v 8.1.1) (GraphPad Software, San Diego, CA, United States) except for Venn diagrams that were created with the jvenn web application⁶ (Bardou et al., 2014). Differences among the four used pipelines, in term of the number of reads assigned, the number of identified phylum and genera as well as their relative abundances were tested using Friedman rank sum test conducted with a Dunn’s correction. Wilcoxon signed rank test was applied when genera relative abundances were compared between two pipelines. Significance level of alpha = 0.05 was used.

A total of 4715000 reads from 40 fecal samples were used to infer and to compare the taxonomical composition identified by the pipeline suggested in QIIME2, Bioconductor, UPARSE/UNOISE and mothur. The processing of the samples took approximately 3 h of computational time in QIIME2, approximately 8 h in Bioconductor, less than 1 h in UPARSE and approximately 9 h in mothur (Table 1). The comparison of the pipelines between Linux and Mac OS showed that QIIME2 and Bioconductor provided identical outcomes, negligible differences in the number of reads assigned at any identity were obtained when using UPARSE or mothur (p > 0.999).

Overall, after filtering and denoising, pipelines in QIIME2 and Bioconductor resulted in a higher number of reads to assign (QIIME2: 3391670; Bioconductor: 3736927; UPARSE: 3173733, regardless the OS; mothur: 3244489, regardless the OS). Among these reads, 84–99% were assigned at the phylum level and 75–99% at the genus level, depending on the pipeline used. Similarly, the proportion of reads assigned to the Bacteria domain, but with no taxonomical assignment at phylum and genus level, differed dramatically among pipelines. Excluding Bioconductor, where the removal of unclassified reads is one step included in the pipeline, the proportion of unclassified reads at the genus level varied from 5% in QIIME2 to 24–29% in UPARSE and mothur (p < 0.001, df = 4, X² = 103). UPARSE-linux provided the same number of unclassified reads at the phylum level as those obtained by using UPARSE-mac (p > 0.999), and the same number of unclassified reads at the genus level as those obtained with mothur, regardless the OS (p > 0.999).

All the used pipelines showed general agreement in the identification of the phyla distribution: the proportion of Firmicutes varied from 46.6 (UPARSE-mac) to 48.3% (QIIME2), Bacteroidetes from 39.7 (QIIME2) to 41.9% (UPARSE-mac), Proteobacteria from 5.7 (QIIME2, Bioconductor, UPARSE) to 6.0% (mothur), Verrucomicrobia from 1.45 (UPARSE, regardless the OS) to 1.90% (Bioconductor), Actinobacteria from 1.44 (UPARSE-linux) to 1.59% (Bioconductor) and, Tenericutes from 0.22 (mothur-mac) to 0.34% (QIIME2) (Figure 2). Cyanobacteria were detected in lower abundance by UPARSE (0.12%, regardless the OS) and mothur (0.18%, regardless the OS) as compared to Bioconductor (1.68%) and QIIME2 (1.21%), whereas Lentisphaerae were identified only by using QIIME2 and Bioconductor (0.47 and 0.50%, respectively). Pairwise comparisons among pipelines revealed widespread differences in the relative abundance of all the phyla (p < 0.0002) (Figure 3). UPARSE (on Linux and Mac OS) or mothur (on Linux and Mac OS) provided comparable values (p > 0.999). Similar results were obtained clustering OTUs at 99% (Supplementary Figure 1).

At the genus taxonomic level, QIIME2 and Bioconductor identified a higher number of bacterial genera than those obtained by UPARSE-linux, UPARSE-mac, mothur-linux, and mothur-mac (p < 0.001, df = 5, X² = 125) (Figure 4A). After singleton removal, 187 genera were identified by QIIME2, 232 by Bioconductor, 120 by UPARSE-linux, 118 by UPARSE-mac, 139 by mothur-linux, and 138 by mothur-mac. Considering all the 316 genera identified, 78 were in common to all the different pipelines, 17 were in common to 5, 19 were in common to 4, 8 were in common to 3, and 87 were in common to 2 out of 6 pipelines. In terms of absolute number of shared genera, QIIME2 and Bioconductor were the two pipelines that provided the most similar results, sharing 147 genera, while UPARSE-mac and QIIME2 were the most different ones, sharing 78 genera. In terms of proportion of shared genera, UPARSE-linux and UPARSE-mac as well as mothur-linux and mothur-mac, were the most producing the overlapping results. Similar results were found after removal of the genera identified in less than five subjects or with less than 10 reads across all samples (Figure 4B).

As the quantitative comparison of all the genera was unfeasible, we arbitrarily selected and compared the 10 most abundant genera identified in each pipeline. All the used pipelines showed general agreement in the identification of the most abundant genera. It has been found that Bacteroides, Faecalibacterium, Alistipes and Subdoligranulum were consistently identified by all the pipelines (Figure 5). Moreover, Blautia and Ruminococcus_1 were concordantly identified by QIIME2, mothur-linux, mothur-mac, UPARSE-linux and UPARSE-mac; CAG-352, Ruminococcaceae ge, Agathobacter and Prevotella 7 were concordantly identified by mothur-linux, mothur-mac, UPARSE-linux and UPARSE-mac and, Parabacteroides, Barnesiella and Ruminococcaceae UCG-002 were concordantly identified by QIIME2 and Bioconductor. Eubacterium coprostanoligenes group was identified by QIIME2 only while Lachnosphira, Ruminococcus 2, and Lachnospiraceae NK4A126 group by Bioconductor only. Pairwise comparisons among pipelines revealed widespread differences in the relative abundance of above-mentioned shared genera. Indeed, pipeline-dependent differences in term of relative abundance was reported for all (p < 0.041) except for Ruminococcus_1 and Ruminococcaceae (p > 0.074). Among the same 10 genera identified as most abundant by UPARSE-linux and UPARSE-mac, a trend forward significance was reported for CAG-352 only (p = 0.056). Similar results were obtained after clustering OTUs at 99% (Supplementary Figure 2).

The 10 genera identified as the most abundant from the different pipelines were selected also when RDP was used as a reference database. Six out of 10 genera were consistently found in all pipelines and the 3 most abundant genera, Bacteroides, Faecalibacterium and Alistipes, have also been identified using SILVA. As in the previous analysis, pairwise comparisons among pipelines revealed widespread differences in the relative abundance of shared genera (p < 0.031) (Supplementary Figure 3).