Clarification of methods for using non-parametric tests to compare groups- vegan package I would like to validate my approach to comparing groups using non-parametric tests. data created below describes the structure of my data. Each var. corresponds to the average concentration of a substance that was measured in an individual (ID) that was randomly sampled from a population of interest (Loc):
set.seed(123)
data <- data.frame(
                   Loc = rep(letters[1:20], each = 20),
                   ID = sample(600,400,replace = F),
                   var1 = rnorm(400),
                   var2 = rnorm(400),
                   var3 = rnorm(400),
                   var4 = rnorm(400),
                   var5 = rnorm(400))

The data frame created above serves only to explain the structure of my data, and allow me to demonstrate how I am doing this analysis (you will not get the same answers I am describing in my questions by using this data, but the process I am using can be replicated with it).
My goal is to determine which locations (Locs) can be differentiated based on the concentrations of these substances, and which substances are contributing to those differences. The substance data does not meet the assumptions required of parametric tests, and my real data set is fairly unbalanced (samples per Loc range from n=6 to n=17). Therefore to answer my questions I have been learning to use the vegan package in R to conduct non-parametric tests, much of which is relatively new to me, which is why I am seeking validation. To perform a PERMANOVA, I started by using the vegdist function to create the distance matrix using Euclidean distance:
vars <- data[,3:7]   
dmat <- vegdist(vars, method = "euclidean")

Next I used the adonis() function to fit the model, and viewed the results via aov.tab:
mod <- Adonis(dmat ~ data$Loc, method = "euclidean") 
mod$aov.tab

The results were significant, indicating that there is a statistically significant difference between at least two Locs in either dispersion or location of their centroids. To rule out over-dispersion I used the betadisper() function to calculate the average distance of each group to the groups centroid, and used both an anova and a permutation test to look for differences:
betad <- betadisper(dmat, data$Loc)
anova(betad)
permutes(betad)

The dissimilarity coefficients did not produce any principle coordinate axes with negative eigenvalues, and both tests gave insignificant results (p<0.05), which suggests that there is indeed support for differences in substance concentrations between the Locs. To further support this idea, I used an anosim on the dissimilarity matrix, which compared the average-ranked dissimilarities between the Locs to the average-ranked dissimilarities within the Locs:
aosim(dmat, data$Loc)

These results came back significant (p>0.05), further suggesting differences between the Locs in terms of substances.
I believe that up to this point I have done everything correctly (correct me if I am wrong, or should have done something differently). What I am really not sure about is the appropriate next steps to take in order to determine where these differences exist. If I were doing a parametric MANOVA, the next thing I would probably do is conduct univariate ANOVAs to compare Locs using each individual var. But in this situation, it is very unclear to me what the appropriate next step would be with a categorical grouping variable and continuous predictors. What can I use to further examine where these differences exist?
 A: You say that your "goal is to determine which locations (Locs) can be differentiated based on the concentrations of these substances, and which substances are contributing to those differences".
Your analysis so far focuses on the first part of your goal - I haven't run through each step, but it looks sensible to me.
Your second goal is to determine which substances contribute to the differences between locations that you have observed. For some context, vegan is normally used for analysis of ecological data, and your question is equivalent to analysis for determining which species contribute most to differences in biological communities.
There are several methods - they generally give similar results so it's up to you what to use. I'll give some names and links:
simper - see https://mb3is.megx.net/gustame/dissimilarity-based-methods/simper
Indicator species analysis - see https://jkzorz.github.io/2019/07/02/Indicator-species-analysis.html
model‐based analysis of multivariate abundance data (mvabund) - https://besjournals.onlinelibrary.wiley.com/doi/10.1111/j.2041-210X.2012.00190.x
