NMDS biplot to show both samples and species using R

Dear Pat,
I have successfully plotted NMDS to show sample sites using following functions.

nmds<-read.table(file=“”, header=T)
plot(nmds$axis1, nmds$axis2)
Then I changed the colors for different sites using functions mentioned in the R tutorial (
Now, I want to show both samples and the species(or OTUs) on the same plot using R functions.

I couldn’t find any related document to do this, so could you please advise me to do this wit R( I know Vegan and Phyloseq packages can be used for this, but I want to plot the NMDS with mothur output values).

Many thanks.

You can certainly do that in R, but I’m afraid it’s a bit beyond the scope of what we have time for in the forum.


as far as I know this isn’t implemented in ggord yet. Here is how I’ve done it

look at this question, it's how I started.
scrs <-, display ='sites'))
scrs <- cbind(scrs, DiagnosisGroup = alpha.clin$DiagnosisGroup)
spec <- envfit(bc.nms, otu, perm=999)

spp.scrs <-, display="vectors")))
spp.scrs <- cbind(spp.scrs, Species =rownames(spp.scrs))
spp.scrs <- cbind(spp.scrs, pVal=spec$vectors$pvals)
# spp.scrs <- cbind(spp.scrs, Taxa=taxa$prefered, phy.col=taxa$phy.col)
spp.scrs <- cbind(spp.scrs, Phyla=taxa$class)

###color vectors by phyla

                 aes(x=0, xend=NMDS1, y=0, yend=NMDS2, color = Phyla),
                 arrow=arrow(length=unit (0.25, "cm")), size=1.5)+
    geom_point(mapping=aes(x=NMDS1, y=NMDS2, shape = as.factor(DiagnosisGroup)), size=2)+
    geom_text(data=scrs, aes(x=NMDS1-.02, y=NMDS2-.02, label=row.names(scrs)), size=4)+