Update: I'm no longer convinced that the ordering is transitive, but it is still described by the pseudomedians. It's just that when you don't have a location shift the pseudomedians need not be additive: you don't necessarily have $m_{XZ}=m_X-m_Z$, so $m_{XZ}=m_{ZY}-m_{YZ}$. I think everything else holds -- because most of the results below are negative rather than positive.
We can think in terms of the estimated pseudo-medians and their standard errors.
If the distributions are continuous and differ only by location shift, not in shape or scale, then the large-sample distributions of the pseudo-medians will be $N(m_i,\sigma^2/n)$ for $m_X$, $m_Y$ and $m_Z$ pseudo-medians and the same $\sigma^2$. Since $\hat m_Z-\hat m_X=(\hat m_Z-\hat m_Y)+(\hat m_Y-\hat m_X)$, the $X:Z$ $p$-value will be smaller than the minimum of the other two if the other two are both less than 0.5 (so that the pseudomedians are ordered the right way).
When the distributions differ in scale or shape it's more complicated, because the $\sigma^2$s will differ. For example, if $Y$ is less variable than $Z$ and $X$ it's a lot easier to get the $p$-values not to be ordered the way you want. For example, consider
r<-replicate(1000,{x<-rnorm(500)
y<-rnorm(500,0.15,s=0.1)
z<-rnorm(500,0.4)
c(wilcox.test(x,y,paired=TRUE)$p.value,
wilcox.test(y,z,paired=TRUE)$p.value,
wilcox.test(x,z,paired=TRUE)$p.value)})
table(r[3,]<pmin(r[1,],r[2,]))
table(r[3,]<pmax(r[1,],r[2,]))
You get quite a few examples where $p_{XY}>p_{XZ}>p_{YZ}$ even with all the $p$-values less than 0.5.
However, if $\hat m_X\sim N(m_X, \sigma^2_X/n)$ in large samples (and similarly for $Y$ and $Z$), we know that the variance of $m_X-m_Y$ is $(\sigma^2_X+\sigma^2_Y)/n$ (and similarly).
Now, $\sigma^2_X+\sigma^2_Z$ is not necessarily less than $\max\{\sigma^2_X+\sigma^2_Y, \sigma^2_Y+\sigma^2_Z\}$, so we might not expect even to have $p_{XZ}<\max \{ p_{XY}, p_{YZ}\}$, but it's going to be a pain to prove.
Simulation to the rescue:
> set.seed(2023-9-18)
> r<-replicate(10000,{x<-rnorm(500)
+ y<-rnorm(500,0.1,s=0.1)
+ z<-rnorm(500,0.2)
+
+ c(wilcox.test(x,y,paired=TRUE, alternative="less")$p.value,
+ wilcox.test(y,z,paired=TRUE, alternative="less")$p.value,
+ wilcox.test(x,z,paired=TRUE, alternative="less")$p.value)})
> which(!(r[3,]<pmax(r[1,],r[2,])))
[1] 1420 1745 2226 3959 5440 6158 9463 9844
Now:
> r[,1420]
[1] 0.4756392 0.5053687 0.6037722
Nope: one is bigger than 0.5
> r[,1745]
[1] 0.6671444 0.6444861 0.6835040
Nope: they are all bigger than 0.5
> r[,2226]
[1] 0.2439124 0.2258590 0.2484954
Success: it is possible to have $p_{XZ}<\max \{ p_{XY}, p_{YZ}\}$. The constraint does seem to be almost true, but not quite.
I think any sharp results will have to be on the Normal-approximation scale rather than on the $p$-value scale.
