The inverse Mellin transform provides another approach that uses residues.
I had previously asked about evaluating a similar series using this method.
For $\operatorname{Re} (s)>2$, we have
$$ \begin{align} \left\{\mathcal{M} \ \frac{1}{\sinh^{2} (\pi x)} \right\}(s) &= \int_{0}^{\infty} \frac{x^{s-1}}{\sinh^{2}(\pi x)} \, dx \\ &= \pi^{-s}\int_{0}^{\infty} \frac{u^{s-1}}{\sinh^{2} u} \, du \\ &= 4 \pi^{-s} \int_{0}^{\infty} u^{s-1} \, \frac{e^{-2u}}{(1-e^{-2u})^2} \, du \\ &= 4 \pi^{-s} \int_{0}^{\infty} u^{s-1} \sum_{n=1}^{\infty} n e^{-2nu} \,du \\ &= 4 \pi^{-s} \sum_{n=1}^{\infty} n \int_{0}^{\infty} u^{s-1} e^{-2nu} \, du \\ &= 4 \pi^{-s} \sum_{n=1}^{\infty}n \, \frac{\Gamma(s)}{(2n)^{s}}\\ &= \pi^{-s} \ 2^{2-s} \ \Gamma(s) \zeta(s-1). \end{align}$$
So according to the Mellin inversion theorem,
$$\frac{1}{\sinh^{2} (\pi x)} = \frac{1}{2 \pi i} \int_{c - i \infty}^{c + i\infty} \pi^{-s} \ 2^{2-s} \ \Gamma(s) \zeta(s-1) x^{-s} \, ds \, , $$ where $c= \frac{5}{2}$.
(Since the Mellin transform of $\frac{1}{\sinh^{2}(\pi z)}$ converges absolutely in the right half-plane $\operatorname{Re}(s) >2$, we could have chosen any value of $c$ greater than $2$.)
Replacing $x$ with $n$ and summing both sides, we get
$$\begin{align} \sum_{n=1}^{\infty} \frac{1}{\sinh^{2} (\pi n)} &= \frac{1}{2 \pi i} \int_{5/2 - i \infty}^{5/2 + i\infty} \pi^{-s} \ 2^{2-s} \ \Gamma(s) \zeta(s-1) \zeta(s) \, ds \\ &= \frac{1}{2 \pi i} \int_{5/2 - i\infty}^{5/2+ i\infty} f(s) \, ds. \end{align}$$
Since the magnitude of $\Gamma(z)$ decays exponentially fast to $0$ as $\Im(z) \to \pm \infty$, we can shift the contour to the line $\operatorname{Re}(s) =1$ and conclude that $$\sum_{n=1}^{\infty} \frac{1}{\sinh^{2} (\pi n)} = \frac{1}{2 \pi i} \left( \operatorname{PV} \int_{-\infty}^{\infty} f(1+it) i \, dt + \pi i \, \text{Res}[f,1] + 2 \pi i \, \operatorname{Res}[f,2] \right),$$ where
$$\begin{align} \operatorname{Res} [f(s),1] &= \lim_{z \to 1} (s-1) \zeta(s) \pi^{-s} 2^{2-s} \Gamma(s) \zeta(s-1) \\ &= (1) \left( \frac{1}{\pi} \Big)(2)(1)\Big(-\frac{1}{2} \right) \\ &= - \frac{1}{\pi } \end{align}$$
and
$$ \begin{align} \operatorname{Res}[f(s),2] &= \lim_{z \to 2} (s-2) \zeta(s-1) \pi^{-s} 2^{2-s} \Gamma(s) \zeta(s) \\ &= (1) \left(\frac{1}{\pi^{2}} \Big) (1)(1) \Big(\frac{\pi^{2}}{6} \right) \\ &= \frac{1}{6}. \end{align}$$
But it can be shown using the functional equation of the Riemann zeta function that $f(s)$ is odd (and purely imaginary) along the line $\operatorname{Re}(s)=1$.
Specifically, $$f(1+it) = 4i \cosh \left(\frac{\pi t}{2} \right)\operatorname{csch}(\pi t) \left|\zeta(it) \right|^{2} \, , \quad t \in \mathbb{R}. $$
Therefore, $$\operatorname{PV} \int_{-\infty}^{\infty} f(1+it) i \, dt =0 $$ and the result follows.
EDIT:
Alternatively, it can be shown that $$\sum_{n=1}^{\infty} \frac{1}{\sinh^{2} (\pi n)} = 4 \sum_{n=1}^{\infty} \frac{n}{e^{2 \pi n}-1}. $$
See this answer. This is the series I had asked about in my question.
Also see this paper for the evaluation of similar sums using the inverse Mellin transform.