The way a computation of an integral along an unbounded contour reduces to a sum of residues is often left unjustified, as is done in the linked MO post. In fact this requires an analysis of the behavior of the integrand around $z=\infty$; for $z^{-2}\tanh^3 z$, one takes a bounded part of the contour (say, with $|z|<R$), makes it closed (say, the boundary of $[-R,R]+i\pi[-N,N]$ slit along the positive real axis, where $N$ is an integer), and shows that all the "extra" things vanish in the limit (as $N,R\to\infty$). In our case this doesn't work because of the blow-up of $e^{-z}$ as $\Re z\to-\infty$. Still this can be fixed (see the answer by @Svyatoslav).
Another approach (cf.) is as follows. For $a,b\geqslant 0$ and $0<\Re s<1$ we have $$\int_0^\infty x^{s-2}(e^{-ax}-e^{-bx})\,dx=\frac{\Gamma(s)}{1-s}(b^{1-s}-a^{1-s})$$ (say, use $e^{-ax}-e^{-bx}=x\int_a^b e^{-xy}\,dy$ and justify the $dy\,dx\mapsto dx\,dy$), then we find $$\int_0^\infty x^{s-2}\tanh x\,dx=\frac{\Gamma(s)}{1-s}\sum_{n=0}^\infty\big(2(4n+2)^{1-s}-(4n)^{1-s}-(4n+4)^{1-s}\big)$$ (use $\tanh x=(1-e^{-2x})^2\sum_{n=0}^\infty e^{-4nx}$ and DCT for termwise integration). The sum equals $2^{2-s}(1-2^{2-s})\zeta(s-1)$ (compute $\sum_{n=1}^\infty$ for $\Re s>2$, and use analytic continuation; the last function of $s$ is entire). Thus we get (still for $0<\Re s<1$) $$\int_0^\infty x^{s-2}(\tanh x-xe^{-x})\,dx=\Gamma(s)\big(f(s)-1\big),\\f(s):=\frac{2^{2-s}}{1-s}(1-2^{2-s})\zeta(s-1).$$
It just remains to take $s\to 0^+$; then $f(s)\to 1$, hence the limit is $$\int_0^\infty\frac{\tanh x-xe^{-x}}{x^2}\,dx=f'(0)=1-12\zeta'(-1)-\frac73\log2.$$
The expression stated in the OP now follows from the relation between $\zeta'(-1)$ and $\zeta'(2)$ one obtains from the functional equation for $\zeta$, as well as the value of $\zeta'(0)$. Note also that $1-12\zeta'(-1)=12\log A$, where $A$ is the Glaisher–Kinkelin constant.