There is no way to do this without extensive own programming. Or there are, but they are either ways you do not want to have, or are not correct. For factor variables the permutation test (anova) only gives the effect of the factor, not its individual levels.
Firstly, permutest function has one type of test that you do not have in anova for constrained ordination: argument by = "onedf" will give you sequential test of one-degree-of-freedom contrasts of factors. However, you probably do not want to have sequential tests that compare each contrast to previous levels, and with default contrasts, the first given effect is the corner-point contrast to the first level of your factor (and the first level will not appear in the results)
Secondly, function coef will give you the regression coefficients of single contrasts for each axis (and what you get will depend on the type of your contrasts), but you should also get the standard deviations of these coefficients to get their significance. Moreover, even if you get the sd's, there is a coefficient for each effect for each axis, but no overall test over all axes. There is function vcov that can get you Canoco-type estimates of standard deviations of coefficients, but the theoretical basis of these sd's is lacking and I would not trust them (but provide them in vegan). See the Example in help(vcov.cca) for calculation of this kind of "t-values".
Third alternative is to use each factor level – Genotype) in your case – as a separate two-level factor variable in model fitting. Easiest way is to use something as data.frame(model.matrix(~ Genotype, data=metadata)) to get the data frame with Genotype split to one-df contrasts. Type of contrasts (see help(contrasts)) will again rule what kind of data frame you get. What ever way you define your contrasts and data frames, you will get one less effects than the levels of your factors. So for five genotypes you can get only four effects. That's a mathematical necessity you cannot circumvent in linear matrix algebra like in constrained ordination.
I suggest you use dbrda instead of capscale. Basically similar, but dbrda handles non-Euclidean distances (such as Cao) whereas capscale is more approximate. Second, the argument for selecting distance function is distance. If you use method like in your example, it will be ignored and Euclidean distances will be used – but that you see in the output that will tell you the type of distances, for instance as Inertia is squared Cao distance, but with method = "cao" it reads Euclidean.
