$\def\des{\operatorname{des}}$Let $\des(\pi)$ be the number of descents of the permutation $\pi$. Then for any permutation $\pi$ in $S_k$, we have
$$\binom{xy+k-\des(\pi)-1}{k} =\sum_{\sigma\tau=\pi}\binom{x+k-\des(\tau)-1}{k}
\binom{y+k-\des(\sigma)-1}{k}.\tag{$*$}\label{star}$$

If $\pi$ is the identity then $\des(\pi)=0$, so the left side of \eqref{star} is $\binom{xy+k-1}{k}$ and on the right, $\tau=\sigma^{-1}$. The OP's identities correspond to the fact that if $k\le 3$ then the number of descents of $\pi$ is the same as the number of descents of $\pi^{-1}$, but this is not true for $k>3$. This also explains the occurrence of the Eulerian numbers as coefficients for $k\le3$. 

As far as I know, the earliest proof of \eqref{star} is in
Bogdan Mielnik and Jerzy Plebański,
*[Combinatorial approach to Baker–Campbell–Hausdorff exponents][1]*,
Annales de l’I. H. P., section A, tome 12, no 3 (1970), pp. 215–254;  further references can be found in Jason Fulman and T. Kyle Petersen, 
*[Card shuffling and P-partitions][2]*, arXiv:2004.01659 [math.CO].



  [1]: http://www.numdam.org/item/?id=AIHPA_1970__12_3_215_0
  [2]: https://doi.org/10.48550/arXiv.2004.01659