I'm creating a violin plot and trying to put brackets comparing the P-Value of MEF group to the other groups (DAY3,DAY6, DAY9), but I have not been able to find a working example that shows how to include multiple P-Values.

Here is the code for the graph:

ggplot(DAY9REoverTimeSimple, aes(x=Day, y=Openness, fill = Day)) +
geom_violin() +
geom_boxplot(width=0.1) +
theme_classic() +
theme(axis.text=element_text(size=14), axis.title=element_text(size=14,face="bold")) +
scale_fill_brewer(palette = "Blues") + 
theme(plot.title = element_text(size = 14, face='bold')) 

Here are the P-Values from the ANOVA I did:

res.aov <- aov(Openness ~ Day, data = DAY9REoverTimeSimple)
ANOVA <- summary(res.aov)
Tukey <- TukeyHSD(res.aov)
> Tukey
Fit: aov(formula = Openness ~ Day, data = DAY9REoverTimeSimple)

$Day
             diff         lwr        upr     p adj
DAY3-MEF    7.0096448  3.16557140 10.8537183 0.0000052
DAY6-MEF    7.6491654  3.76950158 11.5288292 0.0000006
DAY9-MEF    4.8194108  1.09515551  8.5436661 0.0033798