How do you make one factor show as symbol, and another factor as colour in nMDS (vegan)?

Question

I am trying to make an nMDS plot of data with a nested factor. I would like the nMDS to show both factors on one plot by using symbols and colour.

In this reproducible example, if use was nested in moisture, I would like the plot to show Moisture as different symbols, and then Use as different colours.

So far I have figured out this:

library("vegan")
library("BiodiversityR")

data(dune, dune.env)
MDS <- metaMDS(dune, distance="bray", strata=dune.env$Moisture)
MDS

plot(MDS$points[,2], MDS$points[,1], type="n", main="Communities by Use", 
     xlab="NMDS Axis 1", ylab="NMDS Axis 2", xlim=c(-1.5,1.5), ylim=c(-1.5,1.5))
ordisymbol(MDS, dune.env, factor="Use", cex=1.25, rainbow=T, legend=T)

Which gives me the different uses as both different symbols and colours, but shows me nothing about moisture. Is it possible to make it show the different factors instead? I'm assuming it might be somewhere in the MDS$points[,] arguments but I'm not sure what exactly those are doing.

Answer 1

Figured it out by modifying the answer from this question: Plot points of metaMDS

data(dune, dune.env)
dune.MDS <- metaMDS(dune, distance = "bray", strata=dune.env$Moisture)
dune.MDS

pchs<- c(0:5)
gr.moi <- factor(dune.env$Moisture)
gr.use <- factor(dune.env$Use)
col.gr <- c("red", "blue", "purple")


plot(dune.MDS, type = "n", display = "sites")
orditorp(dune.MDS,display="species",col="dark grey",air=0.01)
points(dune.MDS, display = "sites", pch = pchs[gr.moi], col = col.gr[gr.use])
legend("topright", legend=levels(gr.moi), bty = "n", col= c("black"), pch = pchs)
legend("bottomright", legend = levels(gr.use), bty = "n", col = col.gr, pch=c(20),)

And it will produce a lovely plot with symbols and colours exactly how I wanted :)