## 2 Materials and Methods
All samples received were stored and frozen prior to processing. For the hypervariable region, V3-V4 of microbial 16S ribosomal RNA a library extraction was conducted after DNA extraction. Prior to sequencing on an Illumina MiSeq platform (Thermo Fischer Scientific, Waltham, MA, USA), the libraries were tagged with individual index pairs, verified using an Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA), and pooled. The raw reads' quality was evaluated using FastQC and MultiQC [22]. Only sequences with a minimum length of 50 were kept after the low-quality sequences were filtered and trimmed with Trimmomatic [23]. Low-quality calls and the initial 12 base calls were eliminated (Phred score \(<\) 20, sliding window size 5). A Kraken2 was applied to carry out the read categorization [24; 25], with k-mer size 31 against the SSU Ref NR 99 database (release 132) of SILVA [26]. Bracken [27] was used to evaluate the taxa abundances and microbiome composition. The Shannon index was calculated, which is a diversity index used by statistics to summarize the diversity of a population in which each member belongs to a unique group; in our case, the groups are bacteria. Version 4.05 of the R program (R Foundation for Statistical Computing, Vienna, Austria) was applied for statistical analysis. A \(p\)-value of 0.05 or less was considered significant. For 16s rRNA analysis, the linear discriminant analysis (LDA) effect size (LEfSe) method [28] was used, which determines the features most likely to explain differences between the two groups (potential biomarkers for differentiating them). It was performed with default parameters on the Galaxy server of the Huttenhower lab ([[http://huttenhower.sph.harvard.edu/galaxy/](http://huttenhower.sph.harvard.edu/galaxy/)]([http://huttenhower.sph.harvard.edu/galaxy/](http://huttenhower.sph.harvard.edu/galaxy/)), accessed on 3 November 2022). The alpha level for the factorial Kruskal-Wallis test among the groups was 0.05 and the threshold for the logarithmic LDA score for discriminative features was 2. For alpha diversity analysis, the Wilcoxon rank sum test was used between the phylum and the genus level Shannon indexes of the groups. For beta diversity analysis, the principal coordinates analysis (PCoA) based on the Bray-Curtis distance measure was applied, which is a statistical method that converts data on distances between items into map-based visualization of those items. The association between the Shannon index and maximal CCA IMT at two taxonomic levels was analyzed using linear regression including Shannon indexes, age, and gender as covariates. Since most covariates did not show a significant difference in the univariate analysis and also to prevent collinearity, we did not include all covariates in the regression, only age and gender. To find the relationship between the bacterial groups, principal component analyses (PCAs) were performed, which is a statistical technique for reducing the dimensionality of a large dataset containing a high number of dimensions/features per observation and analyze it. The PCA was performed in Matlab programming language (MATLAB R2020b) with the use of the Fathom toolbox ([[https://www.usf.edu/marine-science/research/matlab-resources/fathom-toolbox-for-matlab.aspx](https://www.usf.edu/marine-science/research/matlab-resources/fathom-toolbox-for-matlab.aspx)]([https://www.usf.edu/marine-science/research/matlab-resources/fathom-toolbox-for-matlab.aspx](https://www.usf.edu/marine-science/research/matlab-resources/fathom-toolbox-for-matlab.aspx)), accessed on 3 November 2022). The relative abundances of bacteria between the two groups were compared using permutation tests to compare the medians of the two groups.