Milk Consumption Patterns and Gut Microbiota Composition in Healthy Adults: Insights from a Propensity Score-Matched Cross-Sectional Study

Aim: Milk consumption patterns may be associated with gut microbiota composition, but observational findings can be strongly influenced by demographic and body composition differences among consumers. This observational cross-sectional study evaluated whether questionnaire-derived habitual milk consumption patterns were associated with gut microbiota composition in healthy adults by comparing non-milk consumers, bovine milk consumers, and plant-based beverage consumers. Methods: Baseline cross-sectional data from 265 participants were analyzed: 99 non-milk consumers, 139 bovine milk consumers, and 27 plant-based beverage consumers. Clinical, anthropometric, body composition, and microbiota data were integrated. Alpha diversity was assessed using standard diversity indices, whereas beta diversity was evaluated using Bray–Curtis distances and permutational multivariate analysis of variance. Taxon-level analyses were corrected using the false discovery rate, and pairwise propensity score matching was applied as a sensitivity analysis to evaluate the influence of measured demographic, anthropometric, and body composition imbalance. Results: Crude analyses showed a small beta-diversity difference between bovine milk and plant-based beverage consumers (R2 = 0.0119, p = 0.014), alongside marked differences in age, body mass index, fat mass, and impedance-derived parameters. In matched sensitivity analyses, the crude beta-diversity difference between bovine milk consumers and plant-based beverage consumers was not reproduced in a smaller matched subset (pair-stratified PERMANOVA R2 = 0.0214, p = 0.506). No matched comparison showed significant alpha-diversity differences, and no taxa remained significant after false-discovery-rate correction. Conclusions: Crude beta-diversity differences between bovine milk consumers and plant-based beverage consumers were small and were not reproduced in propensity score-matched sensitivity analyses. These findings should be interpreted cautiously because of categorical exposure assessment, residual imbalance, limited PBBC sample size, dispersion differences in crude analyses, and potential unmeasured confounding.