I'm the maintainer of the ggbiplot pkg, so I dug in and modified the code for ggbiplot() to do what I want here. I added a geom.ind argument that allows geom.ind = c("point", "text"). For testing purposes, this is still on a new geoms branch in the package.

#remotes::install_github("friendly/ggbiplot", ref = "geoms")
# adjust variable names to fold at '_'
vn <- rownames(peng.pca$rotation)
vn <- gsub("_", "\n", vn)
rownames(peng.pca$rotation) <- vn

ggbiplot(peng.pca, 
         choices = 3:4,
         groups = peng$species, 
         ellipse = TRUE, ellipse.alpha = 0.1,
         circle = TRUE,
         var.factor = 4.5,
         geom.ind = c("point", "text"),
         point.size = 2,
         labels = lab, labels.size = 6,
         varname.size = 5,
         clip = "off") +
  theme_minimal(base_size = 14) +
#  theme_penguins("dark") +
  theme(legend.direction = 'horizontal', legend.position = 'top') 

This gives: