As usual check carefully (although there is a code snippet to illustrate).
I seem to have a form of mathematical dyslexia.
Answer: $\left[x^{M}\right]\frac{\Gamma\left(x+a\right)\Gamma\left(x+b\right)}{\Gamma\left(x\right)\Gamma\left(x+a+b\right)}\left(x\right)_{\left(a+b\right)}$
Actually the expression gives the answer for all $M$. Some code is included for that. I can expand it to include M and the binomial summation on a if needed.
The process can probably be specialized to just compute one answer for a particular $M$. Let me know if you need it.
First some standard facts:
Definition 1. The Pochhammer symbol.
$$\begin{align*}\left(x\right)_{n}=x\cdot\left(x+1\right)\ldots\left(x+n-1\right)
\end{align*}$$
Definition 2. Unsigned Stirling Number
$$\left[\begin{array}{c}
n\\
l
\end{array}\right]\equiv\left[x^{l}\right]{\displaystyle \prod_{k=0}^{n-1}\left(x+k\right)=\left[x^{l}\right]x\cdot\left(x+1\right)\ldots\left(x+n-1\right)}=\left[x^{l}\right]\left(x\right)_{n}$$
$$=\left[x^{l}\right]{\displaystyle \sum_{j=0}^{n}}\left[\begin{array}{c}
n\\
j
\end{array}\right]x^{j}$$
We can partition
$$\left(x\right)_{\left(a+b\right)}\rightarrow\left(x\right)_{a}\cdot\left(x+a\right)_{b} \tag{1}\label{1}$$
The formula for all convolutions:
$$(x)_{a}\cdot\left(x\right)_{b} \tag{2}\label{2}$$
The problem can be phrased as the convolution :
$${\displaystyle {\displaystyle \sum_{j=0}^{min\left(a,b,M\right)}}\,\left[\begin{array}{c}
a\\
j
\end{array}\right]\left[\begin{array}{c}
b\\
M-j
\end{array}\right]}={\displaystyle \sum_{j=0}^{min\left(a,b\right)}}\left(\left[x^{j}\right]\left(x\right)_{a}\right)\cdot\left(\left[x^{M-j}\right]\left(x\right)_{b}\right)$$
$$=\left[x^{M}\right]\left(\left(x\right)_{a}\cdot\left(x\right)_{b}\right) \tag{2}\label{3}$$
Remark. The upper limit can be extended but additional terms would be zero.
The intent is to convert (1) to (2). Which can be done using:
http://functions.wolfram.com/06.10.17.0004.02
Substituting our symbols.
$$\left(x\right)_{b}=\frac{\Gamma\left(x+a\right)\Gamma\left(x+b\right)}{\Gamma\left(x\right)\Gamma\left(x+a+b\right)}\left(x+a\right)_{b}$$
Which is so obvious it's laughable.
Thus the answer is:
$$\left[x^{M}\right]\frac{\Gamma\left(x+a\right)\Gamma\left(x+b\right)}{\Gamma\left(x\right)\Gamma\left(x+a+b\right)}\left(x\right)_{\left(a+b\right)}=\left[x^{M}\right]\left(x\right)_{a}\cdot\left(x\right)_{b}$$
Maxima example code for all M:
load ("stirling")$
gamma_expand: true$
gamma(x+5);
p1:pochhammer(x,6);
expand(p1);
p2:pochhammer(x+a,6);
conv(a,b):=pochhammer(x,a+b)*gamma(x+a)*gamma(x+b)/((gamma(x))*gamma(x+a+b));
ratsimp(conv(4,5));
ratsimp(pochhammer(x,4)*pochhammer(x,5));
https://www.researchgate.net/publication/259476881_Convolution_identities_for_Stirling_numbers_of_the_first_kind_via_involution https://www.emis.de//journals/INTEGERS/papers/k9/k9.pdf – rrogers Feb 20 '20 at 23:20