Laurent expansion of the gamma function at negative integers

It is known that the gamma function has the non-positive integers as its simple poles. And we know the residues at these poles.

$$\operatorname*{Res}_{z=-n}\Gamma(z) = \frac{(-1)^n}{n!} .$$

This means, its Laurent expansion around $$z = -n$$ is of the form

$$\Gamma(z) = \frac{(-1)^n}{n!} \frac{1}{z+n} + C_0 + C_1 (z+n) + \cdots .$$

The question is, what is the value of $$C_0$$, or even $$C_1$$?

The expansion of $$\Gamma(z-n)$$ around $$z=0$$ can be obtained from known $$\Gamma(1+z)=1-\gamma z+\gamma_2 z^2+\ldots,\quad\gamma_2=\frac{1}{2}\Big(\gamma^2+\frac{\pi^2}{6}\Big)$$ using properties of $$\Gamma$$ and partial fractions: $$\frac{\Gamma(z-n)}{\Gamma(1+z)}=\prod_{k=0}^{n}\frac{1}{z-k}=\frac{(-1)^n}{n!}\sum_{k=0}^{n}\binom{n}{k}\frac{(-1)^k}{z-k};$$ expanding each term using geometric series, we get $$\frac{\Gamma(z-n)}{\Gamma(1+z)}=\frac{(-1)^n}{n!}\left(\frac{1}{z}+\sum_{r=0}^{\infty}A_r z^r\right),\quad A_r=\sum_{k=1}^{n}\binom{n}{k}\frac{(-1)^{k-1}}{k^{r+1}}$$ (actually $$A_0=\sum_{k=1}^{n}\frac{1}{k}$$ is just $$n$$-th harmonic number; here are further details).
Thus $$C_0=\dfrac{(-1)^n}{n!}(A_0-\gamma)$$ and $$C_1=\dfrac{(-1)^n}{n!}(A_1-\gamma A_0+\gamma_2)$$.