Place p-value at the top of ggplot bar graph using stat_compare_means and facet_wrap, and perform Tukey test on specific comparison

I have two questions related to performing statistics and adding p-values to ggplot graphs.

I know that you can separately generate statistics and add text annotations on ggplot, but it seems easiest to do it all within ggplot if possible. Still, I'm open to other solutions that solve the problem.

I'm using the iris data set. Here is my code.

iris
melt_iris <- melt(iris)

#set specific comparisons
my_comparison_iris <- list(c("setosa", "versicolor"), c("setosa", "virginica"))

#plot data
ggplot(data = melt_iris, aes(x = Species, y = value, fill = Species)) + 
   geom_bar(color = "black", stat = "summary", fun = "mean") +
   geom_point(size = 5) +
   scale_fill_manual(values = c("gray65", "deepskyblue1", "orange")) +
   scale_color_manual(values = c("darkslategray", "darkolivegreen", "darkslateblue", "deeppink4")) +
   theme_bw() +
   theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
         legend.position = "top",
         legend.title = element_text(face = "bold", size = 14),
         legend.text = element_text(size = 12),
         strip.text = element_text(size = 10, face = "bold"),
         legend.background = element_rect(size=.5, linetype="solid", color = "black"),
         legend.box = "vertical") +
   scale_y_continuous(expand = expansion(mult = c(0.05, 0.2))) +
   stat_compare_means(comparisons = my_comparison_iris) +
   stat_compare_means(method = "anova", label.y.npc = "top") + 
   facet_wrap(~ variable, ncol = 4, scales = "free_y") 

First, I cannot get the anova p-value to print on the top of the graph. I am using facet_wrap so I can't give a specific y-value. I tried using label.y.npc = "top" or = 1, which I thought would put the label at the top of each axis. I also added white space using scale_y_continuous and expand, but that didn't help. How can I move the anova p-value to the top of each graph?

Second, I'm unclear on the default method is for stat_compare_means and am having trouble finding a list of accepted options. I think the way I wrote it, it's doing a pairwise t-test or Kruskal-Wallis, which is not corrected for multiple tests. What I would like to do, assuming the overall anova shows significance, is perform multiple pairwise comparisons using the Tukey test. Is that an option?

Thank you for any help you can provide.