RDA visualisation - plot Species abundance according to RDA1

I have sampled 87 plots and identified overall 9 earthworm species. My aim is to investigate changes in species richness and abundance across my site (sampled across 87 plots), which is known to be polluted (mainly by Cu, Zn, and Pb but I also considered 2 composite measure of soil environment extracted by a PCA).

For that I have performed an RDA (species data were hellinger-transformed):

rda(formula = species.hel ~ Cu + Zn + Pb + PCA_axis1 + PCA_axis2, data = env)

the model is significant, I do have an effect of my studied pollutants on species community and the variance explained by my explanatory variables reach 37%. All is ok so far except that now, I would like to plot how individual species abundance varies with RDA1 such as in Oksanen online course (http://cc.oulu.fi/~jarioksa/opetus/metodi/ordination101.html#33):

If my understanding of how RDA proceeds is correct, from the RDA object i should potentially extract (non) linear species response to environmental variable with one model for each species. However, in the summary object of the RDA, I cannot figure out which data to extract. I any of you has a suggestion on that I would be very grateful.

The source code of that lecture is available in github. Just look for the piece of code that produces the graphic that you displayed above, and change the cca() to rda with your model and run the code. The resulting response is expected to be linear, though (your expectation of non-liner response is non-standard).