I am not sure what you have tried, but the "vegan" library (package) has the ability to do rarefaction. So make sure to load library("vegan") and then something like rarecurve(BCI, step = 20, sample = raremax, col = "blue", cex = 0.6) . Documentation for rarefy can be found here. In this tutorial we can see they use rarecurve() for phyloseq data, there they used rarecurve(t(otu_table(ps)), step=50, cex=0.5) but using rarecurve() from vegan should allow you to generate a rarefaction curve from the phyloseq otu_table.