Since the proof is not that difficult let's present the result alluded in the comment above due to Alexandrov (see the monograph by Ross and Shapiro). We present the result only as analytic continuation goes, but it is much stronger and holds for generalized analytic continuation too, and then assuming $S$ satisfies an extra arithmetic condition (which the squares do), one can show that generalized analytic continuation implies analytic continuation in that case, so the OP function doesn't extend beyond the circle even in a weaker sense
Theorem (Alexandrov): Assume $S$ is a set of nonnegative integers st for every $\epsilon>0$ there are $k,m$ positive integers with $k/m < \epsilon$ and with $S$ included in the union of at most $k$ residue classes modulo $m$ except for possibly finitely many elements (that can depend on $m$).
Then if $f(z)=\sum a_nz^n$ is a power series supported on $S$ (ie $a_n =0, n \notin S$) that converges on some disc $D(0,R)$ and if $f$ has analytic continuation at some $|w|=R$ then $f$ has an analytic continuation at all points of modulus $R$ so in particular its radius of convergence is strictly greater than $R$.
Corollary: If $\sum a_nz^n$ supported on $S$ has a radius of convergence precisely $R$, then the circle of radius $R$ is its natural boundary
Lemma: The set of squares satisfies the hypothesis of the theorem.
Proof: for any odd prime $p$ the set of squares is included in $\frac{p+1}{2} \le \frac {2p}{3}$ residues modulo $p$ so for any distinct odd primes $p_1,..p_r$ the set of squares is included in at most $\Pi \frac{p_k+1}{2} \le (\frac{2}{3})^rp_1..p_r$ residues modulo $p_1..p_r$ so the required result for a fixed $\epsilon$ follows for $r$ large enough so $(\frac{2}{3})^r<\epsilon$
Proof of the Theorem: assume $f$ extends at some point $w$, hence it extends on some arc $I$ of the circle containing $w$ of some length $2c>0$; by a rotation $f(\alpha z)$ one can assume wlog that the arc is $0 < \theta < 2c$ (where the circle is parametrized by $\theta \in [-\pi, \pi]$) and then pick $k,m$ as in the hypothesis with $\frac{2\pi k}{m} <c$
Let $q_j, j=1,..k$ the residues mod $k$ that contain $S$ (we can ignore finitely many elements in $S$ since the part of $f$ obtained from them is a polynomial etc, so can assume wlog $S \subset \cup_{b\ge 0, j=1,..k}bm+q_j$) and let $\omega=e^\frac{2\pi i}{m}$ the fundamental root of unity of order $m$.
The crucial fact now is that for each $j=1,..k$ the function $f(\omega^jz)$ extends on the interval $-\frac{2\pi j }{m}< \theta < 2c-\frac{2\pi j }{m}$ which intersects $I$ in at least $0<\theta<c$ since $2c-\frac{2\pi j}{m} \ge 2c-\frac{2\pi k}{m} >c, j=1,..k$, so if we show that we can write $f$ as a linear combination of $f_j(z)=f(\omega^jz), j=1,..k$, it follows that $f$ extends analytically across at least $(-d, 2c)$ where $d=\frac{2\pi }{m}$ and then of course we can iterate this by rotating $f$ so it extends across $(0, 2c+d)$ to get an extension to $(-d, 2c+d)$ etc so in finitely many steps we cover the full circle.
But letting $F_j(z)=\sum_{n=bm+r_j}a_nz^n$ one has that $f(z)=\sum_{j=1}^kF_j(z)$, while $f_j(z)=\sum_{l=1}^j\omega^{jq_l}F_l(z)$
Since trivially the square determinant with terms $\omega^{jq_l}, j=1,..k, l=1,..k$ is non zero as the $\omega^{q_l}$ are distinct for distinct $l$, we can invert the relations above and get $F_l(z)=\sum_{j=1}^k c_{jl}f_j(z)$, hence $f(z)=\sum_{j=1, l=1}^k c_{jl}f_j(z)$ and we are done!
Best Answer
The equation $\Psi(s) = \int_C z \Gamma(z)^2 \zeta(z) s^{-z} dz$ is correct. The inverse Mellin transform equation is about picking up the residues, there's nothing essential about using a straight-line contour. In fact, as you have probably noticed, a straight-line contour won't work for this function. Here's a graph of integrand $z \Gamma(z)^2 \zeta(z) s^{-z}$ for negative $s$
(This was generated by the very helpful complex plotter tool which you can find here)
With this in mind, you can take a contour like the following
I'll go back and check for any computational mistakes later, but you should get something like:
$$\frac{1}{s} + -\frac{1}{2} +\sum_{n=1}^{\infty}\frac{s^{n}}{n!^{2}}\left(\zeta(-n)\left(1+n\ln(s)-2n\psi(1+n)\right)-n\zeta'(-n)\right)$$