Mellin transform of digamma function
This problem is solved quickly if you apply the integral form of digamma function, and then manage to simplify the double improper integral, which require a few tricks. The integral form of digamma that I used to solve the problem was:
$$\psi(t) = \int_0^\infty \left({e^{-s} \over s }- {e^{-s t} \over {1 - e^{-s}}}\right)ds $$
Just combine the last integral form into the definition of Mellin Transform to create a double improper integral:
$$M_x (\psi(x+1)) = \int_0^\infty t^{x-1} \psi(t+1) dt = \int_0^\infty \int_0^\infty t^{x-1} \left({e^{-s} \over s }- {e^{-s t} \over {1 - e^{-s}}}\right) ds dt$$
Combining the expression inside the double integral, and factor the common terms, we obtain:
$$M_x (\psi(x+1)) = \int_0^\infty \int_0^\infty {{t^{x-1} e^{-s} (1 - s e^{-s t} - e^{-s})} \over {s (1 - e^{-s})}} ds dt$$
We begin the tricky part, just as understand that the part of denominator and the numerator can be evaluated as a result of a geometric series, this means:
$$\sum_{n=1}^\infty e^{-n s} = {e^{-s} \over {1 - e^{-s}}}$$
This make the double integral more easy to evaluate, by change part of the function into a power series:
$$M_x(\psi(x+1)) = \sum_{n=1}^\infty \int_0^\infty \int_0^\infty {{t^{x-1} e^{-n s} (1 - s e^{-s t} - e^{-s})} \over {s }} ds dt$$
The integral inside the power series can be factored with two parts, each can be evaluated quickly. Ignoring so far the full expression of the problem to avoid lost of focus, this means the double integral is splitted into:
$$\int_0^\infty \int_0^\infty {{t^{x-1} e^{-n s} (1 - s e^{-s t} - e^{-s})} \over {s }} ds dt =$$ $$ \int_0^\infty \int_0^\infty {{t^{x-1} ( e^{-s n} - e^{-s (n+1)})} \over {s }} ds dt - \int_0^\infty \int_0^\infty {{t^{x-1} e^{-s ( n + t)}} } ds dt $$
The first integral is divergent unless an adequate convergence factor (renormalization) are introduced. Also the first integral can be evaluated according to Frullani's formula:
$${\int_0^\infty {{e^{-a t} - e^{-b t}} \over t }} = {\log \left({b \over a} \right)}$$
Then the first integral should be remarked as the following limit of this renormalization function:
$$\int_0^\infty \int_0^\infty {{t^{x-1} ( e^{-s n} - e^{-s (n+1)})} \over {s }} ds dt = \lim_{\alpha \to 0} \int_0^\infty \int_0^\infty {{t^{x-1} e^{-\alpha (s+t)} ( e^{-s n} - e^{-s (n+1)})} \over {s }} ds dt$$
Just expand the exponentials terms, and the integrals can be evaluated:
$$\int_0^\infty \int_0^\infty {{t^{x-1} ( e^{-s n} - e^{-s (n+1)})} \over {s }} ds dt = \lim_{\alpha \to 0} \int_0^\infty \int_0^\infty {{t^{x-1} e^{-\alpha t} ( e^{-s( n + \alpha)} - e^{-s (n+1+\alpha)})} \over {s }} ds dt = ...$$ $$\rightarrow \lim_{\alpha \to 0} \int_0^\infty {{t^{x-1} e^{- \alpha t} \log \left({{\alpha + n + 1}\over {\alpha + n}}\right)}} dt = \lim_{\alpha \to 0} {\Gamma(x) \over \alpha^x} \log\left({{\alpha + n + 1} \over {\alpha + n}}\right)$$
Since the limit inside the logarithm depends of the series parameter, we should evaluated it first:
$$\lim_{\alpha \to 0} \sum_{n=1}^\infty {{\Gamma (x) \over \alpha^x} \log \left({{\alpha + n + 1} \over {\alpha + n}}\right)} = \lim_{\alpha \to 0} {{\Gamma (x) \over \alpha^x} \log \left({ \prod_{n=1}^\infty {{\alpha + n + 1} \over {\alpha + n}}}\right)}$$
The infinite product inside the logarithm is equal to the Wallis Product:
$$ \lim_{\alpha \to 0} \log \left({ \prod_{n=1}^\infty {{\alpha + n + 1} \over {\alpha + n}}}\right) = \log \left( {\prod_{n=1}^\infty {{n + 1} \over {n}}}\right) = \log \left({\pi \over 2} \right) $$
And the other limit converges if and only if $x<0$:
$$\lim_{\alpha \to 0} {\Gamma (x) \over \alpha^x} \log \left({\pi \over 2}\right) \rightarrow 0 , x<0 $$
This means the first integral is conditionally convergent to zero, if and only if $x<0$.
About the second integral, the integration is straightforward, it just a product of two gamma functions:
$$- \int_0^\infty \int_0^\infty {{t^{x-1} e^{-s ( n + t )}} } ds dt = - n^{x-1} \Gamma(x) \Gamma(1-x) $$
And finally the sum is finite, if we restrict the domain to $1-x>1 \rightarrow x<0$, so:
$$\sum_{n=1}^\infty - n^{x-1} \Gamma(x) \Gamma(1-x) = - \zeta(1-x) \Gamma(x) \Gamma(1-x)$$
Joining alltogether, the Mellin Transform of Digamma Function are:
$$M_x(\psi(x+1)) = - \zeta(1-x) \Gamma(x) \Gamma(1-x)$$
The domain of Mellin Transform of Digamma function can be analytically continued by the functional equation of zeta function, which is defined by: $$ \zeta(s) = 2^{s} \pi^{s-1} \sin \left({\pi s \over 2}\right) \Gamma(1-s) \zeta(1-s) $$
However to justify the result announced by this question, just apply the Euler's Reflection Formula: $$\Gamma(s) \Gamma(1-s) ={ \pi \over \sin\left(\pi s\right)}$$
And finally we get: $$M_x(\psi(x+1)) = - {\pi \over \sin(\pi x)} \zeta(1-x)$$