NGS/operational taxonomic unit (OTUs) data, LC-MS data of silymarin components and of the catabolites were integrated using multiblock sparse partial least-squares discriminant analysis (sPLS-DA), implemented as a DIABLO method in the R package mixOmics (Version 6.26.0). Visualization was provided by methods from the mixOmics package. First, the block.sPLS-DA function was used to determine the optimal number of components based on the performance of the model when considering the centroid distance technique and the lowest balanced error rate, with a fivefold cross-validation (repeated 100 times). The optimized model contains 30 OTUs, three silymarin components for LC-MS and five silymarin microbial catabolites for LC-MS. In addition, the plotDiablo function was used to generate a plot showing the overall correlation between the most discriminant amplicon sequence variants and catabolites, with a network function used to visualize correlations of > 0.6.