#Initial RDA with ALL of the environmental data
mite.spe.rda<-rda(mite.spe.hel~., data=mite.env)
 
#Select significant environmental variables
ordiR2step(rda(mite.spe.hel~1, data=mite.env), 
           scope= formula(mite.spe.rda), direction= "forward", R2scope=TRUE, pstep=1000)
 
#Create a new dataframe with only the significant variables that you identified above
mite.env.signif <- subset(mite.env, 
                          select = c("WatrCont", "Shrub", "Substrate", "Topo", "SubsDens"))
 
#Re-run the RDA with the significant variables and look at the summary
mite.spe.rda.signif=rda(mite.spe~., data=mite.env.signif)
summary(mite.spe.rda.signif, display=NULL)
 
#Find the R2 adjusted of the model with the retained environmental variables 
(R2adj <- RsquareAdj(mite.spe.rda.signif)$adj.r.squared)
 
#Determine the significant canonical (constrained) axes)
anova.cca(mite.spe.rda.signif, step=1000)
anova.cca(mite.spe.rda.signif, step=1000, by="axis")  
 
#Plot the RDA
windows()
plot(mite.spe.rda.signif, scaling=1, main="Triplot RDA - scaling 1", type="none", xlab=c("RDA1"), ylab=c("RDA2"), xlim=c(-1,1), ylim=c(-1,1))
points(scores(mite.spe.rda.signif, display="sites", choices=c(1,2), scaling=1),
       pch=21, col="black", bg="steelblue", cex=1.2)
text(scores(mite.spe.rda.signif, display="species", choices=c(1), scaling=1),
     scores(mite.spe.rda.signif, display="species", choices=c(2), scaling=1),
     labels=rownames(scores(mite.spe.rda.signif, display="species", scaling=1)),
     col="grey", cex=0.8)  
arrows(0,0,
      scores(mite.spe.rda.signif, display="bp", choices=c(1), scaling=1),
      scores(mite.spe.rda.signif, display="bp", choices=c(2), scaling=1),
      col="red")
text(scores(mite.spe.rda.signif, display="bp", choices=c(1), scaling=1)+0.05,
     scores(mite.spe.rda.signif, display="bp", choices=c(2), scaling=1)+0.05,
     labels=rownames(scores(mite.spe.rda.signif, display="bp", choices=c(2), scaling=1)),
     col="red", cex=1)