Based on the answers by @James and @Jyotirmoy Bhattacharya I came up with this solution:
zx <- replicate (5, rnorm(50))
zx_means <- (colMeans(zx, na.rm = TRUE))
boxplot(zx, horizontal = FALSE, outline = FALSE)
points(zx_means, pch = 22, col = "darkgrey", lwd = 7)
(See this post for more details)
If you would like to add points to horizontal box plots, please see this post.