My calculation shows that
$$I(\epsilon)
:= \int_{0}^{\infty} \frac{dx}{\sinh^{2}(\epsilon\sqrt{x^{2}+1})}
= \frac{\pi}{2\epsilon^{2}} - \frac{1}{\epsilon} + \pi \epsilon \sum_{n=1}^{\infty} \frac{1}{(\pi^{2}n^{2} + \epsilon^{2})^{3/2}} \tag{1}. $$
In particular, if we expand the infinite sum on the RHS, we get
$$ I(\epsilon) = \frac{\pi}{2\epsilon^{2}} - \frac{1}{\epsilon} + \sum_{n=0}^{\infty} \binom{-3/2}{n} \frac{\zeta(2n+3)}{\pi^{2n+2}} \epsilon^{2n+1}. $$
Indeed, from the following expansion
$$ \frac{1}{\sinh^{2}z}
= \sum_{n=-\infty}^{\infty} \frac{1}{(z-i\pi n)^{2}}
= \frac{1}{z^{2}} + 2 \sum_{n=1}^{\infty} \frac{z^{2} - \pi^{2}n^{2}}{(z^{2} + \pi^{2} n^{2})^{2}}, $$
it follows from term-wise integration that
\begin{align*}
\int_{0}^{R} \frac{dx}{\sinh^{2}(\epsilon\sqrt{x^{2}+1})}
= \frac{\arctan R}{\epsilon^{2}} + \sum_{n=1}^{\infty}
&\Bigg( \frac{2\epsilon \arctan \left( R\epsilon \big/ \sqrt{\pi^{2}n^{2} + \epsilon^{2}} \right)}{(\pi^{2}n^{2} + \epsilon^{2})^{3/2}} \\
&\quad + \frac{2}{R}\frac{1}{\pi^{2}n^{2}+\epsilon^{2}} \\
&\quad - \frac{2(R^{2}+1)}{R(\pi^{2}n^{2} + (R^{2}+1)\epsilon^{2})} \Bigg). \tag{2}
\end{align*}
Here, term-wise integration is possible from Fubini's theorem together with the following estimate:
$$ \int_{0}^{R} \left| \frac{\epsilon^{2}(x^{2}+1) - \pi^{2}n^{2}}{(\epsilon^{2}(x^{2}+1) + \pi^{2} n^{2})^{2}} \right| = \frac{\arctan \left( R\epsilon \big/ \sqrt{\pi^{2}n^{2} + \epsilon^{2}} \right)}{\epsilon\sqrt{\pi^{2}n^{2}+\epsilon^{2}}}
\lesssim_{\epsilon, R} \frac{1}{n^{2}}. $$
(Notice that the arctan term also contributes to order $n^{-1}$. It means that this argument fails if we consider $R = \infty$. This is why we consider proper integral first.)
Finally, taking $R \to \infty$ to (2) yields (1). When doing this, the only non-trivial calculation is to check that
$$ \lim_{R\to\infty} \sum_{n=1}^{\infty} \frac{2(R^{2}+1)}{R(\pi^{2}n^{2} + (R^{2}+1)\epsilon^{2})} = \frac{1}{\epsilon}. $$
But this follows from the squeezing lemma combined with the following inequality
$$ C(1, R) \leq \sum_{n=1}^{\infty} \frac{2(R^{2}+1)}{R(\pi^{2}n^{2} + (R^{2}+1)\epsilon^{2})} \leq C(0, R), $$
where
\begin{align*}
C(a, R)
&:= \int_{a}^{\infty} \frac{2(R^{2}+1)}{R(\pi^{2}x^{2} + (R^{2}+1)\epsilon^{2})} \, dx \\
&= \frac{\sqrt{R^{2}+1}}{R\epsilon} \frac{\arctan\left( \epsilon\sqrt{R^{2}+1} \big/ a\pi \right)}{\pi/2}.
\end{align*}