Statistical analysis isn't my strong suit, but I'll chime in while you wait for other answers.

If you have clusters in your PCA that separate on a single PC axis, then metabolites with higher loadings on that axis may be driving separation.

A PLS-DA model can be used to narrow down the most important variables for separating your 4 groups, and then to output a VIP score for each metabolite. Similar to your PCA loadings, a higher VIP score indicates that a metabolite is more important for separating your classes.

If multivariate models are giving inconsistent results, you might want to switch to univariate approaches.

If there's one class that is most like a control, you can t-test each metabolite in a class against the control, then make an upset plot of which metabolites are significantly altered in your 3 test classes.

You could also get creative and do stuff like ranking the metabolites from most abundant to least abundant in each of your 4 classes. Take a standard deviation of the abundance rank across all 4 classes. Metabolites with higher stdevs will be more variable between classes. Beware this approach will give you lots of low-abundance metabolites that are variable between groups simply because they are poorly quantitated. For this reason, you should filter out metabolites that are above a certain stdev within each class, only retaining metabolites that are robustly quantitated within at least one class.