I guess what you want is to explore the position/variability of your bivalves and gastropods in relation to multiple environmental variables and each other? Or, do you have only one gradient: grain size particles?

If the first is correct than I would go for a constrained multivariate analysis (either RDA or CCA) depending on how you thin, the environmental variables are related to the distribution of the species.

If the second one is correct, and this piece of code "log(especies1mm[,1:19]+1)" indicates you log(x+1) transformed the abundance of the species and imported it in the PCA. Then, the PCA is trying to find the "direction/variability" between the abundance of the different species. Hence, the plot you look at probably shows the species with highest and lowest abundance in relation to each other along axis1/PC1, this is very similar to an PCoA. I guess this is not what you want?

If you have only a single gradient the most easiest way would be a boxplot or you could regress the abundance of each with the grain size particles.

This latter website https://sites.google.com/site/mb3gustame/, might help you in finding a suitable method to explore your data. Additionally, the vegan package in R is very useful (https://cran.r-project.org/web/packages/vegan/vignettes/intro-vegan.pdf) and the following website gives an overview of some R code (https://www.davidzeleny.net/anadat-r/doku.php/en:ordination).

Good luck.