Is this close to what you need?

The code:

# Load packages
library(tidyverse)
library(ggplot2)

# Make some fake NMDS results:
sampsize <- 10
treatments <- c("start", "end", "1%", "4%")
treatment <- rep(treatments, each = sampsize)
nmds1means <- c(-0.3, 0, 0.3, 0.6)
nmds2means <- c(0.2, -0.7, -0.1, -0.4)
nmds1 <- sapply(nmds1means, FUN = function(x){rnorm(n = sampsize, mean = x, sd = 0.12)}) %>% as.vector
nmds2<- sapply(nmds2means, FUN = function(x){rnorm(n = sampsize, mean = x, sd = 0.2)}) %>% as.vector

df <- data.frame(treatment, nmds1, nmds2)

# Plot as "normal" NMDS plot:
ggplot(data = df, mapping = aes(x = nmds1, y = nmds2)) +
  geom_point()

# Prepare data for summarised NMDS plot:
df2 <- df %>%
  dplyr::group_by(treatment) %>%
  dplyr::summarise(nmds1_mean = mean(nmds1),
                   nmds1_se = sd(nmds1),
                   nmds2_mean = mean(nmds2),
                   nmds2_se = sd(nmds2))

# Plot summarised NMDS plot:
ggplot(data = df2, mapping = aes(x = nmds1_mean, y = nmds2_mean, shape = treatment)) +
  geom_point(size = 4) +
  scale_shape_manual(values = c(15, 16, 17, 18)) +
  geom_errorbar(mapping = aes(ymin = nmds2_mean - nmds2_se, ymax = nmds2_mean + nmds2_se)) +
  geom_errorbar(mapping = aes(xmin = nmds1_mean - nmds1_se, xmax = nmds1_mean + nmds1_se)) +
  lims(x = c(-1, 1), y = c(-1, 1)) +
  xlab("NMDS1") +
  ylab("NMDS2")