$$
f(z) = \frac{\cot z \coth z}{z^{4n-1}} \quad \text{and} \quad \gamma_N \quad \text{the square with corners} \quad (\pm 1 \pm i)\left(N + \frac{1}{2}\right)\pi.
$$
Due to $\cot z = -\frac{2}{z} \sum_{k=0}^{\infty} \frac{\zeta(2k)}{\pi^{2k}} z^{2k}$ and $\coth z = -\frac{2}{z} \sum_{k=0}^{\infty} (-1)^{k} \frac{\zeta(2k)}{\pi^{2k}} z^{2k}$, by the Cauchy product formula,
$$
f(z) = \frac{\cot z \coth z}{z^{4n-1}} = \frac{4}{z^{4n+1}} \sum_{m=0}^{\infty}\sum_{k=0}^{m}(-1)^{k} \frac{\zeta(2k)}{\pi^{2k}} \frac{\zeta(2m-2k)}{\pi^{2m-2k}} z^{2m}.
$$
The term of the form $\text{Coefficient} \cdot \frac{1}{z}$ is obtained for $m=2n$.
Therefore, $\text{res}f\big|_{z=0} = 4 \sum_{k=0}^{2n}(-1)^{k} \frac{\zeta(2k)}{\pi^{2k}} \frac{\zeta(4n-2k)}{\pi^{4n-2k}}$.
Furthermore, $\text{res}f\big|_{z=k\pi} = \lim_{z\to k\pi} \frac{z-k\pi}{\sin z} \frac{\cos z \coth z}{z^{4n-1}} = \frac{\coth k\pi}{(k\pi)^{4n-1}}$ and $\text{res}f\big|_{z=ik\pi} = \lim_{z\to ik\pi} \frac{z-ik\pi}{\sinh z} \frac{\cot z \cosh z}{z^{4n-1}} = \frac{\coth k\pi}{(k\pi)^{4n-1}}$.
Since $\sin z = \sin(x+iy) = \sin x \cosh y + i\cos x \sinh y$, we have $|\sin z|^2 = \sin^2x(\cosh^2y + 1) + \cos^2x \sinh^2y = \sin^2x + \sinh^2y$, which is at least $1$ when $x$ or $y$ is $\pm (N+\frac{1}{2})\pi$.
Considering $\cot^2z = \frac{1}{\sin^2z} - 1$, for all $z\in \gamma_N$, $|\cot z|^2 \leq 1+1$. Also, $z\in \gamma_N \implies iz\in \gamma_N$, so $|\coth z|^2 = |\cot iz|^2 \leq 2$.
Thus, $|f(z)| \leq \frac{2}{|z|^{4n-1}}$ for all $z\in \gamma_N$. This implies $\left|\oint_{\gamma_N}f\,dz\right| \leq 8\left(N+\frac{1}{2}\right)\frac{2}{\left(N+\frac{1}{2}\right)^{4n-1}} \to 0$ as $N\to \infty$.
Therefore, the sum of all residues must be $0$:
$$
4\sum_{k=1}^{\infty}\frac{\coth k\pi}{(k\pi)^{4n-1}} + 4\sum_{k=0}^{2n}(-1)^{k}\frac{\zeta(2k)}{\pi^{2k}}\frac{\zeta(4n-2k)}{\pi^{4n-2k}} = 0
$$
- Proof:
$\frac{\pi}{k}\coth(k\pi) = \sum_{\ell \in \mathbb{Z}}\frac{1}{k^{2}+\ell^{2}} = \frac{1}{k^{2}} + 2\sum_{\ell=1}^{\infty}\frac{1}{k^{2}+\ell^{2}}$
Divide by $(k^{2})^{2n-1}$: $\frac{\pi}{k^{4n-1}}\coth(k\pi) = \frac{1}{k^{4n}} + 2\sum_{\ell=1}^{\infty}\frac{1}{(k^{2})^{2n-1}(k^{2}+\ell^{2})}$
Summing over $k$: $$\pi\sum_{k=1}^{\infty}\frac{\coth(k\pi)}{k^{4n-1}} = \zeta(4n) + 2\sum_{k,\ell=1}^{\infty}\frac{1}{(k^{2})^{2n-1}(k^{2}+\ell^{2})}$$
Due to symmetry, $$\sum_{k,\ell=1}^{\infty}\frac{1}{(k^{2})^{2n-1}(k^{2}+\ell^{2})}$$ can be written as $$\sum_{k,\ell=1}^{\infty}\frac{1}{(\ell^{2})^{2n-1}(k^{2}+\ell^{2})}$$
Hence, $$2\sum_{k,\ell=1}^{\infty}\frac{1}{(k^{2})^{2n-1}(k^{2}+\ell^{2})} = \sum_{k,\ell=1}^{\infty}\frac{1}{(\ell^{2})^{2n-1}(k^{2}+\ell^{2})} + \sum_{k,\ell=1}^{\infty}\frac{1}{(k^{2})^{2n-1}(k^{2}+\ell^{2})}$$
$$= \sum_{k,\ell=1}^{\infty}\frac{1}{(k^{2})^{2n-1}(\ell^{2})^{2n-1}}\frac{(k^{2})^{2n-1}+(\ell^{2})^{2n-1}}{k^{2}+\ell^{2}}$$
where $$\frac{(k^{2})^{2n-1}+(\ell^{2})^{2n-1}}{k^{2}+\ell^{2}} = \sum_{r=1}^{2n-1}(-1)^{r-1}(k^{2})^{2n-1-r}(\ell^{2})^{r-1}$$
So, $$\pi \sum_{k=1}^{\infty}\frac{\coth(k\pi)}{k^{4n-1}} = \zeta(4n) + \sum_{k,\ell=1}^{\infty}\sum_{r=1}^{2n-1}(-1)^{r-1}\frac{1}{(k^{2})^{r}(\ell^{2})^{2n-r}}$$
$$= \zeta(4n) + \sum_{r=1}^{2n-1}(-1)^{r-1}\zeta(2r)\zeta(4n-2r) = \sum_{r=0}^{2n}(-1)^{r-1}\zeta(2r)\zeta(4n-2r).
$$