DNA isolation, library preparation and 16S rRNA gene sequencing
Extracted bacterial DNA was subjected to polymerase chain reaction (PCR) amplification of the V4 region of the 16S rRNA bacterial gene using 2-stage Nextera PCR libraries using the primer pairs 515F (5′- GTG YCA GCM GCC GCG GTA A -3′) and 806R (5′- GGA CTA CNV GGG TWT CTA AT -3′).Raw sample extracts were diluted to 2.5 ng/μL using Tris buffer, and 5 μL were used in first step PCR, together with 5x HOT FIREPol MultiPlex Mix (Solis BioDyne) and 4 μM primer mix (forward and reverse) 515F/806R (Microsynth). First step PCR samples were purified with NGS Clean Beads (Labgene). Bead ratio was 1:1:2, and beads were washed with 75% ethanol, airdried, and resuspended in Tris buffer. In second step PCR, each sample was individually barcoded, using Nextera XT Index Kit v2 (Illumina, San Diego, California) and 5x HOT FIREPol MultiPlex Mix (Solis BioDyne). Second step PCR samples were purified with NGS Clean Beads (Labgene). The final second step PCR products were quantified using a Quant-iT PicoGreen dsDNA Assay Kit (Thermo Fisher Scientific). Amplicons were pooled in equimolar amounts prior to sequencing. The final pool was quantified using a Quant-iT PicoGreen dsDNA Assay Kit (Thermo Fisher Scientific) and Fragment analyzer (Agilent).Subsequent PCR libraries were sequenced on an Illumina MiSeq platform using a v2 500 (2×250 bp read length). Pools were diluted to 9.2 pM and loaded together with 15% PhiX (Illumina, FC-110-3001) to increase the diversity of the run resulting in a raw cluster density of 631 and a cluster passed filter rate of 98%. Produced paired-end reads that passed Illumina’s chastity filter were subjected to demultiplexing and trimming of Illumina adaptor residuals using Illumina’s bcl2fastq software version v2.20.0.422. The quality of the reads was checked with the software FastQC version 0.11.8, and sequencing reads that fell below an average Q-score of 20 or had any uncalled bases (N) were removed from further analysis. The locus-specific V4 primers were trimmed from the sequencing reads with the software cutadapt v3.2. Paired-end reads were discarded if the primer could not be trimmed. Trimmed forward and reverse reads of each paired-end read were merged to reform in silico the sequenced molecule considering a minimum overlap of 15 bases using the software USEARCH version 11.0.667. Merged sequences were again quality filtered, allowing a maximum of one expected erroneous base per merged read. Reads that contain ambiguous bases or were outliers regarding the amplicon size distribution were also discarded. Samples that resulted in less than 5000 merged reads were discarded so as not to distort the statistical analysis. Remaining reads were denoised using the UNOISE algorithm implemented in USEARCH to form amplicon sequencing variants (ASVs), discarding singletons and chimeras in the process. The resulting ASV abundance table was then filtered for possible barcode bleed-in contaminations using the UNCROSS algorithm. ASV sequences were compared to the reference sequences of the RDP 16S database provided by 