Bacterial and fungal PCR, sequencing, and sequence analysis and Taxonomic composition

Bacterial 16 S rRNA genes were amplified using polymerase chain reaction (PCR) targeting the V4 region with dual-barcoded, as per procedure as described in²⁵. Next, amplicons sequenced with an Illumina MiSeq using the 250-bp paired-end kit (v.2). Sequences were denoised, taxonomically classified using Greengenes (v. 13_8) as the reference database, and clustered into 97% similarity operational taxonomic units (OTUs) with the mothur software package (v. 1.39.5) previously described²⁶, following the recommended procedure (https://www.mothur.org/wiki/MiSeq_SOP; accessed August 2018). The resulting dataset had 21257 OTUs (including those occurring once with a count of 1, or singletons). An average of 18383 quality-filtered reads generated per sample. Sequencing quality for R1 and R2 was determined using FastQC 0.11.5.

ITS2 region were sequenced on an Illumina MiSeq (v. 2 chemistry) using the dual barcoding protocol as described²⁵. Primers and PCR conditions used for 16 S rRNA gene and ITS2 sequencing were identical to those previously described²⁷. Bacterial sequences were processed and clustered into operational taxonomic units (OTUs) with the mothur software package (v. 1.39.5)²⁶, following the recommended mothur SOP. Paired-end reads were merged and curated to reduce sequencing error as described in²⁸. The resulting dataset had 3171 OTUs (including those occurring once with a count of 1, or singletons). An average of 9581 quality-filtered reads were generated per sample. Sequencing quality for R1 and R2 was determined using FastQC 0.11.5. Fungal processing pipeline was identical as the one used for bacteria, except for the following differences: (1) paired-end reads were trimmed at the non-overlapping ends, and (2) high quality reads were classified using UNITE (v. 7.1) as described before as the reference database²⁹. A consensus taxonomy for each OTU obtained and the OTU abundances then aggregated into genera. OTU table was rarified to 10000 reads per sample to correct for differences in sequencing depth with rarefy_even_depth function of phyloseq R package³⁰, and alpha diversity indexes (Observed species, Shannon, ACE) were computed from rarified OTU table estimate_richness function of phyloseq R package. The R package vegan was used to compute Beta-diversity matrix from rarified OTU table collapsed at genus level (vegdist function) and to visualize microbiome similaritires with principle coordinate analysis (PCoA) (cmdscale function)³¹. Enterotype classification was performed from the same genus abundance matrix used for PCoA analyses following two different approaches. First, samples were clustered using Jensen-Shannon divergence (JSD) distance and the Partition Around Medoids (PAM) clustering algorithm as described in Aurumugam et al³². Second, samples were clustered from genus abundance data using the Dirichlet Multinomial Mixture (DMM) method of Holmes et al³³. The DMM approach groups samples if their taxon abundances can be modeled by the same Dirichlet-Multinomial (DM) distribution.