There is a good reason to use $z$ instead of $e^{sT}$. Before starting with any analysis, let me remind you that in analysis of signals and systems we are interested in analyzing the frequency spectrum of the signal, i.e. the Laplace transform on the imaginary line $s = jw$. And since your signal $x[n]$ is discrete, then its frequency spectrum its periodic, so its more general define $s=jwT$.
Now, let $X(z) = \mathcal{Z}\{x[n]\}$ of a causal or non-causal discrete signal $x[n]$, i.e.
$$ X(z) = \sum_{n=-\infty}^{\infty} x[n] z^{-n}. $$
Since $z\in\mathbb{C}$ we have $z = |z| e^{j\arg z}$. Without loss of generality we rewrite $|z| = r$ and $\arg z = wT$, i.e. $z=r e^{jwT}$ (note that not necessarily $r=1$). Then
$$ \begin{aligned}
X(z) &= \sum_{n=-\infty}^{\infty} x[n] z^{-n}\\
&= \sum_{n=-\infty}^{\infty} x[n] (r e^{jwT})^{-n}\\
% &= \sum_{n=-\infty}^{\infty} (x[n] r^{-n}) e^{-njwT}\\
&= \sum_{n=-\infty}^{\infty} (x[n] r^{-n}) (e^{jwT})^{-n},
\end{aligned} $$
which implies $X(z) = \left. \mathcal{L}\{x[n]r^{-n}\} \right\rvert_{s=jwT} = \mathcal{F}\{x[n]r^{-n}\}$. As a consequence, $X(z)$ is a Fourier transform more generic than the Fourier transform $X(e^{jwT}) = \mathcal{F}\{x[n]\}$ of our signal of interest.
So, if the convergence radius of $X(z)$ is less than unity then $X(e^{jwT})$ does not exist and therefore its Fourier transform does not either, which represents a problem because there are many signals with this problem of convergence, e.g. non-causal signals such as a digital image filter. Therefore, it is convenient (and even necessary in non-causal signals) to use the Z-transform.
Or informally, use $z$ instead of $\left. e^{sT} \right\rvert_{s=jwT} = e^{jwT}$ whenever you can.
We also recommend to see this link about radius convergence.
At this point, it is clear that the Z-transform has the same objective as the Laplace transform: ensure the convergence of the transform in some region of $\mathbb{C}$, where the Z-transform does it for discrete signals and Laplace transform for continuous signals.
Best Answer
Let $f(t,\mu,\nu)$ be represented by the convolution integral
$$f(t,\mu,\nu)=\int_0^t \tau^{\mu-1}(t-\tau)^{\nu-1}\,d\tau\tag2$$
Note that $f(1,\mu,\nu)=B(\mu,\nu)$, where $B(x,y)$ is the Beta function, which can be represented by the integral
$$B(x,y)=\int_0^1 t^{x-1}(1-t)^{y-1}\,dt$$
Applying the Convolution Theorem of the Laplace Transform to $(2)$ and using $(1)$ reveals that
$$\begin{align} \int_0^\infty f(t,\mu,\nu)e^{-st}\,dt&=\left(\int_0^\infty t^{\mu-1}e^{-st}\,dt\right)\left(\int_0^\infty t^{\mu-1}e^{-st}\,dt\right)\\\\ &=\left(\frac{\Gamma(\mu)}{s^\mu}\right)\left(\frac{\Gamma(\nu)}{s^\nu}\right)\\\\ &=\frac{\Gamma(\mu)\Gamma(\nu)}{s^{\mu+\nu}}\\\\ &=\frac{\Gamma(\mu)\Gamma(\nu)}{\Gamma(\mu+\nu))}\int_0^\infty t^{\mu+\nu-1}e^{-st}\,dt\tag3 \end{align}$$
Taking the inverse Laplace Transform of both sides of $(3)$, we find that
$$f(t,\mu,\nu)=\frac{\Gamma(\mu)\Gamma(\nu)}{\Gamma(\mu+\nu)}\,t^{\mu+\nu-1}\tag4$$
Finally, setting $t=1$ in $(4)$ yields the coveted relationship
as was to be shown!