How to add centroids to an RDA plot

Question

I'd like to replace the arrows on this RDA plot with centroids, something like what's pictured here.

This is the code I currently have which provides me arrows (I guess by default). I have shared our RDA code and I think this is where we might be able to change it from arrows to centroid:

# add arrows for effects of the expanatory variables
arrows(0,0, # start them from (0,0)
       sc_bp[,1], sc_bp[,2], # end them at the score value
       col = "red", 
       lwd = 1,
       length = .1)

(but I share the entire code chunk (below), just in case.

Please note that my data is on fish community (species) and substrate types at 36 sites, I'd like to replace the arrows for substrates with centroids within my RDA.

##Now, the RDA
Y.mat<-Belt_2021_fish_transformed_forPCA #fish community
str(Y.mat)
X.mat<-Reefcheck_2021_forPCA #substrate
str(X.mat)

###Community data has already been transformed with hellinger
##Now, try the RDA
fish_substrate_rda<-rda(Y.mat,X.mat)

```
##Plot
## extract % explained by the first 2 axes
perc_b <- round(100*(summary(fish_substrate_rda)$cont$importance[2, 1:2]), 2)

## extract scores - these are coordinates in the RDA space
sc_si <- scores(fish_substrate_rda, display="sites", choices=c(1,2), scaling=1)
sc_sp <- scores(fish_substrate_rda, display="species", choices=c(1,2), scaling=1)
sc_sp <- sc_sp[c(2,7,8),]
sc_bp <- scores(fish_substrate_rda, display="bp", choices=c(1,2), scaling=1)
sc_bp <- sc_bp[c(2,5,6),]

# Set up a blank plot with scaling, axes, and labels
plot(fish_substrate_rda,
     scaling = 1, # set scaling type 
     type = "none", # this excludes the plotting of any points from the results
     frame = TRUE,
     # set axis limits
     ylim = c(-1.5,0.7),
     xlim = c(-1.5,1.2),
     # label the plot (title, and axes)
     main = "Triplot RDA - scaling 1",
     xlab = paste0("RDA1 (", perc_b[1], "%)"), 
     ylab = paste0("RDA2 (", perc_b[2], "%)") 
)
# add points for site scores
points(sc_si, 
       pch = 21, # set shape (here, circle with a fill colour)
       col = "black", # outline colour
       bg = "steelblue", # fill colour
       cex = 0.7) # size
# add points for species scores
points(sc_sp, 
       pch = 22, # set shape (here, square with a fill colour)
       col = "black",
       bg = "#f2bd33", 
       cex = 0.7)
# add text labels for species abbreviations
text(sc_sp + c(-0.09, -0.09), # adjust text coordinates to avoid overlap with points 
     labels = rownames(sc_sp), 
     col = "grey40", 
     font = 2, # bold
     cex = 0.6)
# add arrows for effects of the expanatory variables
arrows(0,0, # start them from (0,0)
       sc_bp[,1], sc_bp[,2], # end them at the score value
       col = "red", 
       lwd = 1,
       length = .1)
# add text labels for arrows
text(x = sc_bp[,1] -0.01, # adjust text coordinate to avoid overlap with arrow tip
     y = sc_bp[,2] - 0.09, 
     labels = rownames(sc_bp), 
     col = "red", 
     cex = .7, 
     font = 1)
```


I have not found anything online that might help me to accomplish this.