adding points to a surface reverses x and y axis

In trying to replicate some 3D marginal plots of a model, I notice that plot_ly swaps axes when adding points.

When using persp, I get a figure which I consider correct. meshout changes with size, but not distance.

but when mapping size and distance to x and y respectively, axes appear to be switched.

When manually switching them, plot appears fine.

Is this something you guys can look into? Below is a fully working example.

library(plotly)
library(ggplot2)
library(reshape2)
library(rgl)

N <- 50
B0 <-  1.2
B1 <-  1.5
B2 <- -0.5

set.seed(357)
count <- rpois(N, 10)

distance <- seq(-0.5, 0.5,,length(count))
size <- (log(count) - B0 - B2 * distance) / B1
xy <- data.frame(count, size, distance)

mdl <- glm(count ~ size * distance, family = poisson(link = log), data = xy)
summary(mdl)

ggplot(xy, aes(x = size, y = count, color = distance)) +
  geom_point()

npts <- 50
s.size <- seq(from = min(xy$size), to = max(xy$size), length.out = npts)
s.dist <- seq(from = min(xy$distance), to = max(xy$distance), length.out = npts)
mesh <- expand.grid(
  size = s.size,
  distance = s.dist
)
mesh$count <- predict(mdl, newdata = mesh, type = "response")
head(mesh)

meshout <- as.matrix(
  dcast(mesh, formula = size ~ distance, value.var = "count")[, -1]
)

pmat <- persp(x = s.size, y = s.dist, z = meshout, # size and distance mapped OK
      theta = 0, phi = 30,
      col = "grey90", border = "grey50")
points(trans3d(x = xy$size, y = xy$dist, z = xy$count, pmat = pmat))

plot3d(x = xy$size, y = xy$distance, z = xy$count) # check, distance and size are plotted correctly
# https://stackoverflow.com/questions/5468280/scale-a-series-between-two-points#comment49916403_5468527
range02 <- function(x, newMin, newMax){ (x - min(x))/(max(x)-min(x)) * (newMax - newMin) + newMin }

xy$scaled_size <- range02(xy$size, newMin = 0, newMax = 50)
xy$scaled_distance <- range02(xy$distance, newMin = 0, newMax = 50)

# size and distance swapped
plot_ly(z = ~ meshout) %>% 
  add_surface(name = "Incorrect mapping of x and y") %>%
  add_markers(data = xy, x = ~scaled_size, y = ~scaled_distance, z = ~count,
              marker = list(size = 5, color = "black"), name = "Incorrect mapping of x and y") 

# when switching mapping of x and y, plots fine
plot_ly(z = ~ meshout) %>% 
  add_surface(name = "Incorrect mapping of x and y") %>%
  add_markers(data = xy, x = ~scaled_distance, y = ~scaled_size, z = ~count,
              marker = list(size = 5, color = "black"), name = "Incorrect mapping of x and y")

Start with the reproducible R example using plot_ly(), add_surface(), and add_markers(), and compare its axis mapping with persp() and plot3d(). Trace how the surface coordinates are established relative to marker x and y values; done means the surface and added points use the same size and distance axes without manually swapping them.