How to use Euler-Maclaurin summation to show the following relationship over Gamma function

I think you can get this from the Bernoulli integral for the $\Gamma$ function: $$ \Gamma(z) = \int_0^\infty x^{z-1} e^{-x} dx. $$

If we use the first Euler Maclaurin summation and swap the integral and summation, we get: $$ \int_a^b f(x) dx = \sum_{a\leq k\leq b} f(k) - \frac{f(a) + f(b)}{2} - \sum_{1\leq i\leq m} \frac{B_{2i}}{(2i)!} f^{(2i-1)}(N) - R(m). $$

We then truncate the last expansion term summation leaving the simple (0-th order "expansion"): $$ \int_a^b f(x) dx \approx \sum_{a\leq k\leq b} f(k) - \frac{f(a) + f(b)}{2}. $$

For our $\Gamma$ function integral, we let $f(x)=x^{z-1}e^{-x}$ to get: $$ \Gamma(z) = \int_0^\infty x^{z-1} e^{-x} dx \approx \sum_{0\leq k\leq \infty} k^{z-1} e^{-k} - \lim_{a\downarrow 0}\frac{a^{z-1}}{2e^a} - \lim_{b\to\infty} \frac{b^{z-1}}{2e^b} = \sum_{0\leq k\leq \infty} k^{z-1} e^{-k}. $$

Note that for the above to hold, we require $z>0$ so that $f(x)$ is bounded above or $z$ such that $f$ is dominated by a function $g$ with a finite integral $\int_0^\infty g(x) dx$.

Now rewrite $\binom{n+\alpha}{n}$: $$ \binom{n+\alpha}{n} = \frac{\Gamma(n+\alpha+1)}{\Gamma(n+1)\Gamma(\alpha+1)} \approx \frac{\sum_{0\leq k\leq \infty} k^{n+\alpha} e^{-k}}{\Gamma(\alpha+1)\sum_{0\leq k\leq \infty} k^{n} e^{-k}} = \frac{k^\alpha}{\Gamma(\alpha+1)} \frac{\sum_{0\leq k\leq \infty} k^{n} e^{-k}}{\sum_{0\leq k\leq \infty} k^{n} e^{-k}} = \frac{k^\alpha}{\Gamma(\alpha+1)}. $$