How to automatically adjust the width of each facet for facet_wrap?

You can adjust facet widths after converting the ggplot object to a grob:

# create ggplot object (no need to manipulate boxplot width here. 
# we'll adjust the facet width directly later)
p <- ggplot(Data,
       aes(x = trait, y = mean)) +
  geom_boxplot(aes(fill = Ref,
                   lower = mean - sd, 
                   upper = mean + sd, 
                   middle = mean, 
                   ymin = min, 
                   ymax = max),
               lwd = 0.5,
               stat = "identity") +
  facet_wrap(~ SP, scales = "free", nrow = 1) +
  scale_x_discrete(expand = c(0, 0.5)) + # change additive expansion from default 0.6 to 0.5
  theme_bw()

# convert ggplot object to grob object
gp <- ggplotGrob(p)

# optional: take a look at the grob object's layout
gtable::gtable_show_layout(gp)

# get gtable columns corresponding to the facets (5 & 9, in this case)
facet.columns <- gp$layout$l[grepl("panel", gp$layout$name)]

# get the number of unique x-axis values per facet (1 & 3, in this case)
x.var <- sapply(ggplot_build(p)$layout$panel_scales_x,
                function(l) length(l$range$range))

# change the relative widths of the facet columns based on
# how many unique x-axis values are in each facet
gp$widths[facet.columns] <- gp$widths[facet.columns] * x.var

# plot result
grid::grid.draw(gp)

While u/z-lin's answer works, there is a far simpler solution. Switch from facet_wrap(...) to use facet_grid(...). With facet_grid, you don't need to specify rows and columns. You are still able to specify scales= (which allows automatic adjustment of axis scales for each facet if wanted), but you can also specify space=, which does the same thing, but with the scaling of the overall facet width. This is what you want. Your function call is now something like this:

ggplot(Data, aes(x = trait, y = mean)) +
    geom_boxplot(aes(
        fill = Ref, lower = mean-sd, upper = mean+sd, middle = mean, 
        ymin = min, ymax = max),
        lwd = 0.5, stat = "identity") +
    facet_grid(. ~ SP, scales = "free", space='free') +
    scale_x_discrete(expand = c(0, 0.5)) +
    theme_bw()

Some more description of layout of facets can be found here.

As @cdtip mentioned, this does not allow for independent y scales for each facet, which is what the OP asked for initially. Luckily, there is also a simple solution for this, which utilizes facet_row() from the ggforce package:

library(ggforce)

# same as above without facet_grid call..
p <- ggplot(Data, aes(x = trait, y = mean)) +
  geom_boxplot(aes(
    fill = Ref, lower = mean-sd, upper = mean+sd, middle = mean, 
    ymin = min, ymax = max),
    lwd = 0.5, stat = "identity") +
  scale_x_discrete(expand = c(0, 0.5)) +
  theme_bw()

p + ggforce::facet_row(vars(SP), scales = 'free', space = 'free')