16S rRNA gene sequencing and sequence data processing

The V3 and V4 hypervariable regions of the 16S rRNA gene were amplified by PCR using barcoded universal primers (Supplementary Table ⁵). Sequencing was performed using the Illumina MiSeq platform (Illumina, Inc. San Diego, CA, USA), according to the manufacturer’s instructions. Raw sequence data were analyzed using the QIIME2 (v.2020.11) bioinformatics pipeline. The 300 bp paired-end reads obtained from the Illumina MiSeq platform for each sample were demultiplexed to attribute sequence reads to the appropriate samples and joined. The sequence reads were denoised and dereplicated into amplicon sequence variants (ASVs) using the DADA2 tool, which also filters chimeras. Each read sequence was trimmed to 388 bp. A total of 20,220,280 sequences of the 16S rRNA gene and 7,100 features were generated from 268 swab samples, with a mean frequency of 75,448 sequences per sample. A feature table, equivalent to the ASV table generated using QIIME2, was generated for all samples with a mean frequency of 2,847. The feature table was used for taxonomic classification, alpha and beta diversity analyses, and differential abundance measurements in the different experimental groups. Taxonomy was assigned to each ASV using the SILVA (version 138) database and a fitted classifier classification-sklearn method. Species-level assignments were performed using the BLASTn software (https://blast.ncbi.nlm.nih.gov/). The highest percentage of identity and expectation values were considered when selecting significant BLAST hits.