Adding a bit of diversity to this discussion we suppose that we seek to prove that
$$\coth x = \sum_{n=-\infty}^\infty \frac{x}{x^2+\pi^2 n^2}.$$
Rather than use the standard technique (absolutely nothing wrong with it) of applying a circular contour to
$$ f(z) = \pi \cot(\pi z) \frac{x}{x^2+\pi^2 z^2}$$
we use Mellin transforms, which are admittedly a bit more complicated in this case (than summing the residues at the two poles of $f(z)$ on the imaginary axis at $\pm ix/\pi$), but perhaps have some complex variable didactic value.
Rewrite this as
$$\coth x = \frac{1}{x} + 2 \sum_{n=1}^\infty \frac{x}{x^2+\pi^2n^2}
= \frac{1}{x} + 2x \sum_{n=1}^\infty \frac{1}{x^2+\pi^2n^2}
= \frac{1}{x} + 2x \sum_{n=1}^\infty \frac{1}{n^2} \frac{1}{(x/n)^2+\pi^2}.$$
Now the inner sum, call it $S(x),$ is harmonic and may be evaluated by inverting its Mellin transform.
Recall the Mellin transform identity for harmonic sums with base function $g(x)$, which is
$$\mathfrak{M}\left(\sum_{k\ge 1}\lambda_k g(\mu_k x); s\right) =
\left(\sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s}\right) g^*(s)$$
where $g^*(s)$ is the Mellin transform of $g(x).$
In the present case we have
$$\lambda_k = \frac{1}{k^2}, \quad
\mu_k = \frac{1}{k} \quad \text{and} \quad
g(x) = \frac{1}{x^2+\pi^2}.$$
The Mellin transform of $g(x)$ is
$$\int_0^\infty \frac{1}{x^2+\pi^2} x^{s-1} dx,$$
which we evaluate using a semicircular contour in the upper half plane (no convergence issues as is easily verified), getting
$$g^*(s) (1 + e^{\pi i(s-1)})= 2\pi i
\operatorname{Res}\left(\frac{1}{z^2+\pi^2} z^{s-1}; z = +i\pi\right).$$
Using the branch of the logarithm with arguments from zero to $2\pi$, this becomes
$$g^*(s) = \frac{2\pi i}{1- e^{\pi i s}}\frac{1}{2\pi i} (i\pi)^{s-1}
= \pi^{s-1} \frac{e^{i(\pi/2)(s-1)}}{1- e^{i\pi s}}
= - \pi^{s-1} \frac{i e^{i(\pi/2)s}}{1- e^{i\pi s}}\\
= - \pi^{s-1} \frac{i}{e^{-i(\pi/2)s}- e^{i(\pi/2)s}}
= \frac{1}{2} \pi^{s-1} \frac{2i}{e^{i(\pi/2)s}- e^{-i(\pi/2)s}}
= \frac{1}{2} \frac{\pi^{s-1}}{\sin((\pi/2) s)}. $$
Now in the present case we have
$$\sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s} =
\sum_{k\ge 1} \frac{1}{k^2} \frac{1}{k^{-s}} = \zeta(2-s).$$
It follows that the Mellin transform $Q(s)$ of $S(x)$ is given by
$$ Q(s) = \frac{1}{2} \frac{\pi^{s-1}}{\sin((\pi/2) s)} \zeta(2-s).$$
The Mellin inversion integral for an expansion about zero is
$$\int_{1/2-i\infty}^{1/2+i\infty} Q(s)/x^s ds,$$
shifted to the left. The poles from the sine term are at the even integers, with those in the right half plane being canceled by the trivial zeros of the zeta function.
We have
$$\sum_{q\ge 0} \operatorname{Res}(Q(s)/x^s; s=-2q)
= \sum_{q\ge 0} \frac{(-1)^q}{\pi^{2q+2}} \zeta(2+2q) x^{2q}\\
= \sum_{q\ge 0} \frac{(-1)^q}{\pi^{2q+2}}
\frac{(-1)^q B_{2q+2} (2\pi)^{2q+2}}{2(2q+2)!} x^{2q}
= \sum_{q\ge 0} \frac{2^{2q+2} B_{2q+2}}{2(2q+2)!} x^{2q}.$$
We thus have for the initial sum
$$\frac{1}{x} + 2x S(x) =
\frac{1}{x} + 2x \sum_{q\ge 0} \frac{2^{2q+2} B_{2q+2}}{2(2q+2)!} x^{2q}
= \frac{1}{x} + \frac{1}{x} \sum_{q\ge 0} \frac{2^{2q+2} B_{2q+2}}{(2q+2)!} x^{2q+2}.$$
Recall that the generating function of the Bernoulli numbers is precisely
$$\frac{x}{e^x-1}.$$
Adapting the expression from the harmonic sum calculation to match this generating function (note that the odd index $B_m$ are zero when $m> 1$) we get
$$\frac{1}{x} - \frac{1}{x} + 1 + \frac{1}{x} \sum_{m\ge 0} \frac{B_m}{m!} (2x)^m.$$
Simplify one last time, getting
$$ 1 + \frac{1}{x}\frac{2x}{e^{2x}-1}
= \frac{e^{2x}+1}{e^{2x}-1}
= \frac{e^x+e^{-x}}{e^x-e^{-x}} = \coth x,$$
QED.
Just to recapitulate, the direction of this computation was as follows. We started by conjecturing that a certain harmonic sum represents the hyperbolic cotangent in a neighborhood of zero. Then we computed the Mellin transform of this sum. We inverted that transform for an expansion about zero. The expansion that we obtained was precisely the generating function of the Bernoulli numbers with two exceptions. Substituting that generating function into the sum and correcting for the exceptions we obtained the defining equation of the hyperbolic cotangent in terms of exponentials.
