How to add influential factors in PCA with Bray-Curtis Index without using the BiodiversityR package in R?

I am trying to add loading, or the influential species for the PCOA plot that I have. I calculated the Bray Curtis distance matrix to plot the oordination. Now, I want to add the loadings to determine the most influential species. However, I am unable to use the BiodiversityR package as Mac OS requires to install the tcl-tk. I have been unable to figure out how after trying for three days. So I'm wondering if there is an alternate package or method to figure out the loading to add to the plot. I have to code for what needs to be done and commented how to proceed if the BiodiversityR package had worked.

library(vegan)
set.seed(111)
sp1 <- rnorm(72, mean = 4, 1)
sp2 <- rnorm(72, mean = 2, 1)
sp3 <- rnorm(72, mean = 3, 1)
sp4 <- rnorm(72, mean = 9, 1)

sp.abd <- data.frame(sp1, sp2, sp3, sp4)


species.db <- vegdist(sp.abd, method = "bray")
species.db <- vegdist(sp.abd, method = "bray", upper = TRUE, diag = TRUE)
species.db[is.na(sp.abd)] <- 0
species.pcoa <- cmdscale(species.db, eig = TRUE, k = 3)

speciesREL <- sp.abd
for(i in 1:nrow(sp.abd)){
  speciesREL[i, ] = sp.abd[i, ] / sum(sp.abd[i, ])
}

library(BiodiversityR)
species.pcoa <- add.spec.scores(species.pcoa, speciesREL, method = "pcoa.scores")


##PLOT PCOA
explainvar1 <- round(species.pcoa $eig[1] / sum(species.pcoa $eig), 3) * 100
explainvar2 <- round(species.pcoa $eig[2] / sum(species.pcoa $eig), 3) * 100
explainvar3 <- round(species.pcoa $eig[3] / sum(species.pcoa $eig), 3) * 100
sum.eig <- sum(explainvar1, explainvar2, explainvar3)

df1 <- data.frame(species.pcoa$points)

rda.plot <- ggplot(df1, aes(x=X1, y=X2)) + 
  geom_point(aes(size = 3, alpha = 0.5)) +
  geom_hline(yintercept=0, linetype="dotted") +
  geom_vline(xintercept=0, linetype="dotted") +
  coord_fixed()  +
  theme_classic()

rda.plot

##Adding lines
#df2 <- data.frame(species.pcoa$cproj)

# rda.biplot <- rda.plot +
#   geom_segment(data=df2, aes(x=0, xend=X1, y=0, yend=X2), 
#                color="black", arrow=arrow(length=unit(0.01,"npc"))) +
#   geom_text(data=df2, 
#             aes(x=X1,y=X2,label=rownames(df2),
#                 hjust=0.5*(1-sign(X1)),vjust=0.5*(1-sign(X2))), 
#             color="black", size=4) + 
# rda.biplot 

For example, without any transformation adding the "most influential" species looks something like this: