*2.4. Targeted 16S rRNA Gene Sequencing*

The analytical methods of 16S rDNA analysis were established in a previous study [18]. By referencing Illumina's recommended protocols (https://support.illumina.com/downl oads/16s\_metagenomic\_sequencing\_library\_preparation.html; accessed on 10 December 2020), we performed library construction and amplification of the 16S rRNA gene. In summary, we used the forward and reverse primers 341F and 805R with Illumina overhang adapter sequencing to amplify the V3–V4 region of the bacterial 16S rRNA gene. A Nextera XT Index kit was then used to adjust the dual-index barcodes to the targets in the amplicon and the Illumina sequencing adapters. The quantity and quality of data in the sequenced library were assessed using a QSep100 analyzer (BiOptic, Taipei, Taiwan). Moreover, by using a MiSeq Reagent kit v3, high-throughput sequencing was performed on an Illumina MiSeq 2000 sequencer.

The bioinformatics analytical process was conducted following the workflow described by Callahan et al. [19]. First, by using the R package DADA2 (v 1.14.1), the filtered reads were managed.Taxonomy assignment was administered using the SILVA database (v128) with a minimum bootstrap confidence level of 80 [20]. Multiple sequence alignment of the structural variants was processed with DECIPHER (v2.14.0), and a phylogenetic tree was built from the alignment using phangorn (v2.5.5) [21]. The count table, taxonomy assignment results, and phylogenetic tree were consolidated into a phyloseq object, and community analyses were created by phyloseq (v1.30.0) [22]. One-way ANOVA followed by the Bonferroni post hoc test were utilized to handle multiple comparison analysis. The

analytical process of alpha-diversity and beta-diversity are listed below. The phyloseq package was used to calculate the alpha-diversity. For beta-diversity, principal coordinate analysis (PCoA) was performed on UniFrac distances, and the adonis and betadisper functions from the vegan package (v2.5.6) were used to analyze the dissimilarity of composition among high- and low-consumption groups. The groups were compared with α = 0.05 (Kruskal–Wallis and Wilcoxon tests). The UniFrac package (v1.1) was used to compare the community dissimilarity between groups, demonstrated as UniFrac distances [23]. GraPhlAn [24] helped us to perform the enrichment analysis between the groups, which were analyzed using the linear discriminant analysis (LDA), effect size, (LEfSe) method, and a logarithmic LDA score of more than 2 [25] and were then visualized as a cladogram.
