Best option for 3D scatterplot with different plot symbols

I would like to make a 3D scatterplot with different shape symbols to represent different categories of data. In 2D, this is straightforward using the pch argument in the plot function, and I would like to extend it to the rgl plot3d function

However, as discussed in this thread,

rgl: plot3d with "extended" plotting symbols

rgl::plot3d returns a single symbol regardless of pch argument, while pch3d doesn't embed the symbols in 3d axes. The best option suggested in the thread uses text3d and generates a 3D plot with no axes labels, e.g. for all_symbol a vector of n1 0's, n2 1's, n3 2's etc

library(rgl)
rgl::open3d()
for(i in 1:49){
    rgl::text3d(scores.df$PC1[i], scores.df$PC2[i], scores.df$PC3[i], text = intToUtf8(all_symbol[i]), cex = 2, usePlotmath = TRUE)
}
rgl::box3d()

generates no axes labels or tick marks.

Is there some graphics library that will give me the same basic 3d plotting functionality of rgl while allowing me to do the rather straightforward tasks such as using symbols for classes of data and labeling axes?

Don't use box3d(), use decorate3d() if you want axes with labels. You didn't include any data in your question, so I'll illustrate with fake data:

library(rgl)
fake <- matrix(rnorm(3*49), ncol = 3)
open3d()
pch3d(fake, pch = 1:49, cex = 0.3)
decorate3d()