You already have an answer based on your particular approach, so let me indicate a slightly different and more general approach. The answer is in 3 parts, first is about how differentiation under integral sign works in the holomorphic setting, next is the application to the Gamma function, and finally generalizing to the Mellin transform
Holomorphic Differentiation under Integral:
The nice thing about complex analysis is that by Cauchy's integral formula, if you have a local bound on a holomorphic function, then you also have bounds on all of its derivatives. As a result, the hypotheses for differentiation under the integral sign can be weakened slightly:
Let $(X,\mathfrak{M},\mu)$ be a measure space, $U\subset \Bbb{C}$ open, and $f:X\times U\to \Bbb{C}$ a function such that for each $z\in U$, we have $f(\cdot, z)\in L^1(\mu)$, and for each $x\in X$, $f(x,\cdot)$ is a holomorphic function on $U$. If for each $a\in U$ there is a disc $\Delta_a$ around $a$ and contained in $U$, and there is an integrable $g_a:X\to [0,\infty]$ such that for all $x\in X,z\in \Delta_a$, we have $|f(x,z)|\leq g(x)$, then the function $F:U\to \Bbb{C}$ defined as
\begin{align}
F(z):=\int_Xf(\cdot,z)\,d\mu\equiv \int_Xf(x,z)\,d\mu(x)
\end{align}
is holomorphic in $U$, and for every $n\geq 0$, we have
\begin{align}
F^{(n)}(z)&=\int_X\frac{\partial^n f}{\partial z^n}(x,z)\,d\mu(x).
\end{align}
Proof 1: Using DCT, Cauchy’s inequalities, and MVT:
Fix a point $a\in U$, a closed disk $\overline{D}_{2r}(a)$ around $a$ and contained in $\Delta_a$ (from the hypothesis). Then,
for all $x\in X$ and $z\in D_r(a)$, we have $\frac{\partial f}{\partial z}(x,z)=\frac{1}{2\pi i}\int_{|w-a|=2r}\frac{f(x,w)}{(w-z)^2}\,dw$, by Cauchy’s integral formula. So we have the bound $\left|\frac{\partial f}{\partial z}(x,z)\right|\leq \frac{1}{2\pi}\cdot \frac{g_a(x)}{r^2}\cdot 2\pi (2r)=\frac{2g_a(x)}{r}$. This is the proof for Cauchy’s inequality.
For all $z\in D_r(a)$, we have
\begin{align}
\left|\frac{F(z)-F(a)}{z-a}-\int_X\frac{\partial f}{\partial z}(x,a)\,d\mu(x)\right|&\leq\int_X\left|\frac{f(x,z)-f(x,a)}{z-a}-\frac{\partial f}{\partial z}(x,a)\right|\,d\mu(x).\tag{$*$}
\end{align}
The integrand on the right converges pointwise to $0$, as $z\to a$, simply by definition of the derivative. On the other hand, by the triangle inequality, mean-value inequality, and the above Cauchy inequality, we have
\begin{align}
\left|\frac{f(x,z)-f(x,a)}{z-a}-\frac{\partial f}{\partial z}(x,a)\right|&\leq
\left|\frac{f(x,z)-f(x,a)}{z-a}\right|+
\left|\frac{\partial f}{\partial z}(x,a)\right|\\
&\leq \sup_{\zeta\in D_r(a)}\left|\frac{\partial f}{\partial z}(x,\zeta)\right|+ \left|\frac{\partial f}{\partial z}(x,a)\right|\\
&\leq \frac{2g_a(x)}{r}+\frac{2g_a(x)}{r}\\
&=\frac{4g_a(x)}{r}.
\end{align}
So, we have found an $L^1$ dominating function, so we can apply DCT to $(*)$. This completes the proof for $n=1$.
By arguing similarly (induction if you want to be super formal), you can conclude the same for all $n\geq 0$.
Proof 2: Using Fubini
Alternatively we can argue as follows: since each $f(x,\cdot)$ is holomorphic, we can write for small circular loops $\gamma$ contained in the open set $U$ and $z$ lying inside the circular loop $\gamma$,
\begin{align}
F(z)&=\int_X\left[\frac{1}{2\pi i}\int_{\gamma}\frac{f(x,w)}{w-z}\,dw\right]\,d\mu(x)\\
&=\frac{1}{2\pi i}\int_{\gamma}\int_X\frac{f(x,w)}{w-z}\,d\mu(x)\,dw\\
&=\frac{1}{2\pi i}\int_{\gamma}\frac{F(w)}{w-z}\,dw
\end{align}
where the second equality follows from Fubini's theorem, which we can use because of the estimate $|f(x,w)|\leq g(x)$ (we make the extra assumption that $f$ is measurable on the product and $\mu$ is $\sigma$-finite so Fubini can be applied). This shows $F$ is indeed holomorphic, and it follows immediately (by putting $(w-z)^{n+1}$ in the denominator and running the argument in reverse) that its derivatives are as claimed.
Application to Gamma Function:
For any $x>0$, the function $t\mapsto t^{x-1}e^{-t}$ is in $L^1((0,\infty))$, because near the origin, it behaves like $t^{x-1}$, which is integrable near the origin by an elementary argument; on the other hand near $\infty$, the exponential term is going to make everything decay nicely (in fact, you should convince yourself that $\int_0^{\infty}t^ne^{-t}\,dt = n! <\infty$ for every $n\geq 0$).
In other words, for any $a\in \Bbb{C}$ with positive real part, we can find $0<\delta<\text{Re}(a)$. As a result, in this region, we can take $g(t)=t^{\delta-1}e^{-t}$ as the dominating function. This shows the Gamma function is holomorphic in this region (and since $a$ was arbitrary, $\Gamma$ is holomorphic in the whole right half plane), and also that differentiation under the integral sign is valid.
Generalization to Mellin Transform:
Given a Locally Lebesgue integrable function $f:(0,\infty)\to \Bbb{C}$, we define the Mellin-transform of $f$ by the formula
\begin{align}
(Mf)(s):=\int_0^{\infty}x^{s-1}f(x)\,dx,
\end{align}
and this is of course defined only for those $s\in \Bbb{C}$ such that the function $x\mapsto x^{s-1}f(x)$ is (Lebesgue) integrable on $(0,\infty)$. For convenience, let us denote this set of points by $U_f$. Here are some simple observations:
If a point $s_0=\sigma_0+it_0$ lies in the domain $U_f$, then the whole vertical line $\sigma=s_0+it$, $t\in \Bbb{R}$ lies inside $U_f$ (this is simply because $|x^{s-1}f(x)|=
x^{\text{Re}(s)-1}|f(x)|$).
Let $\alpha,\beta\in\Bbb{R}$, and suppose $f(x)=O(x^{-\alpha})$ as $x\to 0^+$ and $f(x)=O(x^{-\beta})$ as $x\to \infty$. Then, $x^{\sigma-1}f(x)= O(x^{\sigma-1-\alpha})$ as $x\to 0^+$ and $x^{\sigma-1}f(x)=O(x^{\sigma-1-\beta})$ as $x\to \infty$. So, provided $\sigma-1-\alpha>-1$ and $\sigma-1-\beta<-1$ (or equivalently, $\alpha<\sigma<\beta$), it follows that $x\mapsto x^{\sigma-1}f(x)$ is Lebesgue integrable on $(0,\infty)$. In other words, assuming $f$ obeys some reasonable decay conditions, the entire vertical strip $S_{\alpha,\beta}:=\{\alpha<\text{Re}(s)<\beta\}$ lies inside the domain $U_f$.
The Mellin transform of a function $f$ which satisfies the decay conditions above is holomorphic in the entire strip $S_{\alpha,\beta}$ (why? what is the dominating function which allows us to invoke the theorem above)?
Now, observe that if we let $f(x)=e^{-x}$, then its Mellin transform is precisely the Gamma function. Also, in this case we have $\alpha=0,\beta=\infty$, which is why the Gamma function can be defined by this integral formula on the whole right half plane and is holomorphic there.
