Currently, vegan::betadisper is used for the ordination, via principal coordinates. The figure is hand-made via scores and approximate KDE-estimated probability contours.

The default plot method for betadisper is interesting, and likely simpler to interpret when there are many groups / lots of overlap. The ellipse can be adjusted to any approximate probability contour.

par(mar=c(1,1,3,1))
plot(d.betadisper, hull=FALSE, ellipse=TRUE, col=cols)

vegan::metaMDS provides an alternative ordination, with additional functionality to overlay environmental variables or rotation of axes.

m <- metaMDS(d.sub[, d.mu.wide.vars], distance = 'gower')

par(mar=c(1,1,3,1))
o <- ordiplot(m, type='n', axes=FALSE)
points(o, 'sites', cex=1, col=cols[as.numeric(d.sub$.id)], pch=16)

ordisurf(m, d.sub$`Effective Precipitation (mm)`, add=TRUE, col='black')

# location / labeling of centroids must be done manually (?)

box()