I found a solution thanks mainly to @Teodor (and a bit of ChatGPT). My code is quite similar but still somewhat different. Also, there is a lot of customization concerning my plots following the solution. As I oriented myself very closely to a graph from a scientific publication, I decided to post it in case someone wants the code for a (hopefully) publication-worthy looking plot (Fig. 1). I also attached the code for a simple version of the same graph but using ordibar to generate error bars that take covariance into account as proposed by @JariOksanen (Fig. 2).

Teodor's approach:

library(vegan)
library(dplyr)
library(ggplot2)
library(ggvegan)
library(ggrepel)
library(ggpubr)
#Load data
data <- read.csv2("comSFdata.csv", header = TRUE, check.names = FALSE)
#Create new dataframe w/o 1. column
com = data[,2:ncol(data)]
#Turn data from dataframe to matrix format
mcom = as.matrix(com)
#Run NMDS - set.seed for reproducibility
set.seed(123)
nmds1 <- metaMDS(mcom, distance = "bray", k = 2, trymax = 10000)

#Prepare data for, and calculate means + se --> Summarize into single points
nmds_coords <- as.data.frame(scores(nmds1, display = "sites"))
nmds_coords['Treatment'] <- data['Treatment']
nmds_summary <- nmds_coords %>%
  group_by(Treatment) %>%
  summarise(
    NMDS1_mean = mean(NMDS1),
    NMDS2_mean = mean(NMDS2),
    NMDS1_sd = sd(NMDS1),
    NMDS2_sd = sd(NMDS2),
    NMDS1_se = NMDS1_sd / sqrt(n()),
    NMDS2_se = NMDS2_sd / sqrt(n()))

#Prepare data for 2. plot
fort <- fortify(nmds1)

#Generate 2 plots:
p1 <- ggplot(data = nmds_summary, mapping = aes(x = NMDS1_mean, y = NMDS2_mean, colour = Treatment, shape = Treatment)) +
#Customizes samples & plot error bars
  geom_errorbar(mapping = aes(xmin = NMDS1_mean - NMDS1_se, xmax = NMDS1_mean + NMDS1_se), width = 0.05, colour = "black", size = 1) +
  geom_errorbar(mapping = aes(ymin = NMDS2_mean - NMDS2_se, ymax = NMDS2_mean + NMDS2_se), width = 0.14383561643835616438356164383562, colour = "black", size = 1) +
  geom_point(size = 6, alpha = 1) +
  scale_colour_manual(values = c("#C77CFF", "#00BFC4", "#F8766D", "#7CAE00")) +
  scale_shape_manual(values = c(16, 16, 15, 17)) +
#Set axis range and force rendering of objects "outside" the range
  coord_cartesian(ylim = c(-0.4, 0.4), xlim = c(-1.15, 1.15), clip = 'off') +
#Remove legend
  theme(legend.position="none") +
#Add custom legend
  annotate("text", x = 1.1, y = 0.365, size = 4, fontface = "bold", hjust = 1, colour = "#C77CFF",
           label = ("Control start")) +
  annotate("text", x = 1.1, y = 0.34, size = 4, fontface = "bold", hjust = 1, colour = "#00BFC4",
           label = ("Control end")) +
  annotate("text", x = 1.1, y = 0.315, size = 4, fontface = "bold", hjust = 1, colour = "#F8766D",
           label = ("1% substrate addition")) +
  annotate("text", x = 1.1, y = 0.29, size = 4, fontface = "bold", hjust = 1, colour = "#7CAE00",
           label = ("4% substrate addition")) +
#Generates coordinate system axis at the origin
  geom_abline(intercept = 0, slope = 0, linetype = "dotted", linewidth = 0.65, colour = "black") +
  geom_vline(aes(xintercept = 0), linetype = "dotted", linewidth = 0.7125, colour = "black") +
#Removes grids & background & adds a frame
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        panel.border = element_rect(colour = "black", fill = "NA", linewidth = 1, linetype = "solid")) +
#Add title as label so it is in the foreground
  annotate("label", x = -1.15, y = 0.36, size = 5.5, label.size = 0, fontface = "bold", hjust = 0, colour = "black",
           label = ("(a) Slump Floor (SF)")) +
#Add custom axis titles
  xlab("NMDS1") + ylab("NMDS2") +
#Customize left plot axis
  theme(axis.ticks.length = unit(-3, "mm"),
        axis.ticks = element_line(linewidth = 1.05),
        axis.title.x = element_text(color="black", size=14, face="bold", vjust = -0.6),
        axis.title.y.left = element_text(color="black", size=14, face="bold", vjust = 10),
#Set custom margin for better scale readability           
        plot.margin = margin(10, 10, 10, 35)) +
  scale_x_continuous(expand = c(0, 0), breaks = c(-1, -0.5, 0, 0.5, 1), labels = c("", "", "","",""), sec.axis = sec_axis(~ ., breaks = c(-1, -0.5, 0, 0.5, 1), labels = c("", "", "","",""))) +
  scale_y_continuous(expand = c(0, 0), breaks = c(-0.2, 0 , 0.2), labels = c("", "",""), sec.axis = sec_axis(~ ., breaks = c(-0.2, 0 , 0.2), labels = c("", "",""))) +
#Set custom x-axis labels
  annotate("text", x = -0.95, y = -0.425, size = 4.5, fontface = "bold", hjust = 1, colour = "black",
           label = ("-1")) +
  annotate("text", x = -0.43, y = -0.425, size = 4.5, fontface = "bold", hjust = 1, colour = "black",
           label = ("-0.5")) +
  annotate("text", x = 0.025, y = -0.425, size = 4.5, fontface = "bold", hjust = 1, colour = "black",
           label = ("0")) +
  annotate("text", x = 0.57, y = -0.425, size = 4.5, fontface = "bold", hjust = 1, colour = "black",
           label = ("0.5")) +
  annotate("text", x = 1.025, y = -0.425, size = 4.5, fontface = "bold", hjust = 1, colour = "black",
           label = ("1")) +
#Set custom y-axis labels  
  annotate("text", x = -1.2, y = 0.4, size = 4.5, fontface = "bold", hjust = 1, colour = "black",
           label = ("0.4")) +
  annotate("text", x = -1.2, y = 0.205, size = 4.5, fontface = "bold", hjust = 1, colour = "black",
           label = ("0.2")) +
  annotate("text", x = -1.2, y = 0, size = 4.5, fontface = "bold", hjust = 1, colour = "black",
           label = ("0")) +
  annotate("text", x = -1.2, y = -0.195, size = 4.5, fontface = "bold", hjust = 1, colour = "black",
           label = ("-0.2")) +
  annotate("text", x = -1.2, y = -0.395, size = 4.5, fontface = "bold", hjust = 1, colour = "black",
           label = ("-0.4"))

p2 <- ggplot() +
#Add custom x-axis label
  geom_point(data = subset(fort, score == 'species'),
             mapping = aes(x = NMDS1, y = NMDS2), alpha = 0) +
#Set axis range and force rendering of objects "outside" the range
  coord_cartesian(ylim = c(-0.4, 0.4), xlim = c(-1.15, 1.15), clip = 'off') +
  geom_segment(data = subset(fort, score == 'species'),
#Adds & customizes vectors for PLFAs
               mapping = aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2),
               arrow = arrow(length = unit(0, "npc"), type = "closed"),
               colour = "black", size = 0.8, alpha = 1)+
#Enable custom label colours and set to bold text
  geom_text_repel(data = subset(fort, score == 'species'),
                  mapping = aes(label = label, colour = label, x = NMDS1 * 1.025, y = NMDS2 * 1.025),
                  alpha = 1, size = 3, fontface = "bold") +
#Remove legend
  theme(legend.position="none") +
#Add custom legend
  annotate("text", x = 1.1, y = 0.365, size = 4, fontface = "bold", hjust = 1, colour = "black",
           label = ("Non-specific")) +
  annotate("text", x = 1.1, y = 0.34, size = 4, fontface = "bold", hjust = 1, colour = "#ff0000",
           label = ("Gram-positive bacteria")) +
  annotate("text", x = 1.1, y = 0.315, size = 4, fontface = "bold", hjust = 1, colour = "#ff8800",
           label = ("Gram-negative bacteria")) +
  annotate("text", x = 1.1, y = 0.29, size = 4, fontface = "bold", hjust = 1, colour = "#00b300",
           label = ("Fungi")) +
#Customize label colours
  scale_colour_manual(values=c("i15:0" = "#ff0000", "a15:0" = "#ff0000", "i16:0" = "#ff0000", "16:1ω7c" = "#ff8800", "16:1ω7t" = "#ff8800", "16:1ω5c" = "#ff8800", "16:0" = "black", "10Me 16:0" = "#ff0000", "i17:0" = "#ff0000", "a17:0" = "#ff0000", "17:1ω7c" = "#ff8800", "cy17:0" = "#ff8800", "17:0" = "black", "10Me 17:0" = "#ff0000", "18:2ω6c" = "#00b300", "18:1ω7c" = "#ff8800",  "18:1ω7t" = "#ff8800", "18:0" = "black", "10Me 18:0" = "#ff0000",  "cy19:0" = "#ff8800", "20:0" = "black")) +
#Generates coordinate system axis at the origin
  geom_abline(intercept = 0, slope = 0, linetype = "dotted", linewidth = 0.65, colour = "black") +
  geom_vline(aes(xintercept = 0), linetype = "dotted", linewidth = 0.7125, colour = "black") +
#Removes grids & background & adds a frame
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        panel.border = element_rect(colour = "black", fill = "NA", linewidth = 1, linetype = "solid")) +
#Add title below grid customization so it is in the foreground
  annotate("label", x = -1.15, y = 0.36, size = 5.5, label.size = 0, fontface = "bold", hjust = 0, colour = "black",
           label = ("(b) Slump Floor (SF)")) +
#Customize right plot axis
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.length = unit(-3, "mm"),
        axis.ticks = element_line(linewidth = 1.05),
        axis.title.x = element_text(color="black", size=14, face="bold", vjust = -0.6),
#Set custom margin for better scale readability        
        plot.margin = margin(10, 10, 10, 10)) +
  scale_x_continuous(expand = c(0, 0), breaks = c(-1, -0.5, 0, 0.5, 1), labels = c("", "", "","",""), sec.axis = sec_axis(~ ., breaks = c(-1, -0.5, 0, 0.5, 1), labels = c("", "", "","",""))) +
  scale_y_continuous(expand = c(0, 0), breaks = c(-0.2, 0 , 0.2), labels = c("", "",""), sec.axis = sec_axis(~ ., breaks = c(-0.2, 0 , 0.2), labels = c("", "",""))) +
#Set custom x-axis labels
  annotate("text", x = -0.95, y = -0.425, size = 4.5, fontface = "bold", hjust = 1, colour = "black",
           label = ("-1")) +
  annotate("text", x = -0.43, y = -0.425, size = 4.5, fontface = "bold", hjust = 1, colour = "black",
           label = ("-0.5")) +
  annotate("text", x = 0.025, y = -0.425, size = 4.5, fontface = "bold", hjust = 1, colour = "black",
           label = ("0")) +
  annotate("text", x = 0.57, y = -0.425, size = 4.5, fontface = "bold", hjust = 1, colour = "black",
           label = ("0.5")) +
  annotate("text", x = 1.025, y = -0.425, size = 4.5, fontface = "bold", hjust = 1, colour = "black",
           label = ("1"))

#Generate multi-panel plot with 2 columns & add/customize labels
ggarrange(p1, p2, ncol = 2, nrow = 1, widths = c(1.06, 1), align = "h")

JariO's approach:

library(vegan)
library(dplyr)
library(ggplot2)
library(ggvegan)
library(ggrepel)
library(ggpubr)
#Load data
data <- read.csv2("comSFdata.csv", header = TRUE, check.names = FALSE)
#Create new dataframe w/o 1. column
com = data[,2:ncol(data)]
#Turn data from dataframe to matrix format
mcom = as.matrix(com)
#Run NMDS - set.seed for reproducibility
set.seed(123)
nmds1 <- metaMDS(mcom, distance = "bray", k = 2, trymax = 10000)
#Prepare data for plotting
fort <- fortify(nmds1)
#Add workaround for adding back Treatment groups
fort['Treatment'] <- wa['Treatment']
#Remove rows without Treatment info
newfort <- na.omit(fort)

label <- c("#C77CFF","#C77CFF","#C77CFF","#C77CFF","#00BFC4","#00BFC4","#00BFC4","#00BFC4","#F8766D","#F8766D","#F8766D","#F8766D","#7CAE00","#7CAE00","#7CAE00","#7CAE00")
label2 <- c("#C77CFF","#00BFC4","#F8766D","#7CAE00")

#Ordibar
ordiplot(nmds1, display = "sites", type = "none", main = "Ordibar NMDS")
ordibar(nmds1, display = "sites", groups = data$Treatment, col = label2, length = 0.2)
points(newfort[, 3], newfort[, 4], pch = 16, cex = 1.5, col = label)
ordiellipse(nmds1, groups = data$Treatment, draw = "lines", col = label2, alpha = 1)