MultcompLetters4 based on Tukey test gave the wrong letters to my plot

Question

I did a tukey test and with multcompLetters4 I implemented the connecting letter display to my boxplot. But the letters are not correct, looking at the means of the groups.

I calculated the means in R and i had big differences between groups, when they had the same letter. I also had groups with different letters but means that are quit similar. This is my code:

# Load libraries
library(ggplot2)
library(readxl)
library(ggplot2)
library(readxl)
library(dplyr)
library(dunn.test)
library(multcomp)
library(multcompView)
rm(list = ls())


# Read data
data <- read_excel("Data_for_boxplots.xlsx", sheet = "Nannochloropsis")
data <- na.omit(data)

# Set the order of levels for 'application method'
data$`application method` <- factor(data$`application method`, levels = c("control", "S 0.01", "S     0.05", "S 0.1", "S 0.5", "I 0.01", "I 0.05", "I 0.1", "I 0.5", "F 0.01", "F 0.05", "F 0.1", "F 0.5"))

#anova
model <- aov (data$`lesion diameter`~ data$`application method`, data= data)
summary (model)

#tukey
tukey <- TukeyHSD(model, conf.level=.95)
print(tukey)

#cld
cld <- multcompLetters4(model,tukey)
print(cld)


#table with factors and 3rd quantile
tk <- group_by(data, `application method`)%>%
  summarise(mean=mean(data$`lesion diameter`), quat = quantile(data$`lesion diameter`, probs = 0.75)) %>%
  arrange(desc(mean))


#extracting the compact letter display and adding to the tk table
cld <- as.data.frame.list(cld$`data$\`application method\``)
tk$cld <- cld$Letters
print(tk)

# Plot the data with only x and y axes visible
p <- ggplot(data, aes(x = `application method`, y = `lesion diameter`, fill = `application method`)) +
  geom_boxplot(fill = c("#A6A6A6", "#6baed6", "#3182bd", "#08519c", "#08306b", "#a1d99b", "#74c476", "#31a354", "#006d2c", "#fcbba1", "#fc9272", "#fb6a4a", "#de2d26")) +
  geom_point(position = position_jitter(width = 0.25, height = 0.01), size = 1.5) +  # Add jitter to the y-position
  geom_text(data = tk, aes(label = cld, x = `application method`, y = quat),
        vjust = -12.5, hjust = 0.5, size = 8)+
  labs(y = "Lesion diameter [cm]", x = "") + 
  scale_y_continuous(breaks = c(0,  0.2,  0.4,  0.6,  0.8,  1,1.2,1.4,1.6,1.8,2,2.2,2.4,2.6,2.8,3),
                 labels = c("0",  "0.2", "0.4",  "0.6",  "0.8",  "1",  "1.2","1.4","1.6","1.8","2","2.2","2.4","2.6","2.8","3")) +
  theme(axis.line = element_line(colour = "black"),# Set color of axes lines
    axis.title.y = element_text(size=17.5, face = "bold"),
    axis.text.y = element_text(size=12.5, face = "bold"),
    panel.background = element_blank(), #remove panel background
    panel.border = element_blank(),  # Remove panel borders
    panel.grid.major = element_blank(),  # Remove major gridlines
    panel.grid.minor = element_blank(),  # Remove minor gridlines
    axis.text.x = element_text(size=12.5, face = "bold", angle = 0, hjust = 0.5)) +  # Rotate x-axis labels$
  guides(fill = FALSE)

# Print the plot
print(p)