Let $f(x) = \tanh(x)\tanh(2x)$. Rewrite the integral at hand as
$$I
= \int_0^\infty \frac{f(x)}{x^2}dx
= \frac12 \int_{-\infty}^{\infty} \frac{f(x)}{x^2} dx
= \frac12 \lim_{k\to\infty}
\int_{-k\pi}^{k\pi} \frac{f(x)}{x^2} dx\\
$$
For each positive integer $k$, consider the rectangular region
$$D_k = \big\{ u + iv \in \mathbb{C} : |u| \le k\pi, 0 \le v \le k\pi \big\}$$
and the contour integral over its boundary:
$$C_k \stackrel{def}{=} \int_{\partial D_k} \frac{f(z)}{z^2} dz
= \left(
\int_{-k\pi}^{k\pi}
+ \int_{k\pi}^{k\pi(1+i)}
+ \int_{k\pi(1+i)}^{k\pi(-1+i)}
+ \int_{k\pi(-1+i)}^{-k\pi}
\right) \frac{f(z)}{z^2} dz
$$
$C_k$ split into 4 pieces, one for each edges of the $D_k$. It is not hard to show
as $k \to \infty$, the contribution from the three edges (the top, left and bottom) goes to zero. This implies
$$I = \frac12 \lim_{k\to\infty} \int_{\partial D_k} \frac{f(z)}{z^2} dz\tag{*1}$$
Let $\phi = e^{2z}$, we have
$$\begin{align}f(z)
&= \tanh(z) \tanh(2z) = \left(\frac{\phi - 1}{\phi+1}\right)\left(\frac{\phi^2-1}{\phi^2+1}\right) = \frac{(\phi-1)^2}{\phi^2+1} = 1 - \frac{2\phi}{\phi^2+1} \\
&= 1 - \frac{1}{\cosh 2z}\end{align}$$
Recall the well known? expansion of $\cot z$
$$\cot z = \sum_{n=-\infty}^\infty \frac{1}{z - n\pi}$$
We have
$$\begin{array}{rrl}
& \frac{1}{\sin z}
&= \frac{1+\cos z}{\sin z} - \frac{\cos z}{\sin z} = \cot\frac{z}{2} - \cot z = \sum_{n=-\infty}^\infty \frac{(-1)^n}{z - n\pi}\\
\implies
& \frac{1}{\cos z}
&= -\frac{1}{\sin(z - \frac{\pi}{2})}
= \sum_{n=-\infty}^\infty \frac{(-1)^{n-1}}{z - (n+\frac12)\pi}\\
\implies
& \frac{1}{\cosh z}
&= \frac{1}{\cos(-i z)} =
-i \sum_{n=-\infty}^\infty \frac{(-1)^n}{z - (n+\frac12)\pi i}\\
\implies
& f(z)
& = 1 + \frac{i}{2} \sum_{n=-\infty}^\infty \frac{(-1)^n}{z - (n+\frac12)\frac{\pi}{2} i}\tag{*2}
\end{array}$$
Inside $D_k$, $f(z)$ has poles at $(n + \frac12)\frac{\pi}{2} i$ for each $n \in \mathbb{N}$. We can evaluate the contour integral in $(*1)$ by summing the residues at these poles. As a result,
$$I = \frac12 (2\pi i)(\frac{i}{2})\sum_{n=0}^\infty \frac{(-1)^n}{((n+\frac12)\frac{\pi}{2}i)^2} = \frac{8K}{\pi}$$
where $\;\displaystyle K = \sum_{n=0}^\infty \frac{(-1)^n}{(2n+1)^2}\;$
is the Catalan's constant.
Update
If one want to minimize the use of complex analysis, an alternate starting point is
the expansion of $f(z)$ in $(*2)$. Let $\alpha_n = (n+\frac12)\frac{\pi}{2}$ and notice
$\alpha_{-1-n} = -\alpha_n$, we have
$$f(x)
= 1 + \frac{i}{2}\sum_{n=0}^\infty (-1)^n\left(\frac{1}{x -\alpha_n i} - \frac{1}{x + \alpha_n i}\right)
= 1 - \sum_{n=0}^\infty \frac{(-1)^n \alpha_n}{x^2 + \alpha_n^2}
$$
Since $f(0) = 0$, we have
$$\frac{f(x)}{x^2}
= \frac{f(x)-f(0)}{x^2}
= \frac{1}{x^2} \sum_{n=0}^\infty (-1)^n \alpha_n\left( \frac{1}{\alpha_n^2} - \frac{1}{x^2 + \alpha_n^2}\right)
= \sum_{n=0}^\infty\frac{(-1)^n}{\alpha_n(x^2+\alpha_n^2)}
$$
If you look at the terms of this series in units of pair, we find the value of any pair
is always non-negative. This means we can integrate this series pair by pair and get:
$$\int_0^\infty \frac{f(x)}{x^2}dx = \sum_{n=0}^\infty \frac{(-1)^n}{\alpha_n}\int_0^\infty \frac{dx}{x^2+\alpha_n^2}
= \frac{\pi}{2}\sum_{n=0}^\infty \frac{(-1)^n}{\alpha_n^2}
= \frac{8}{\pi}\sum_{n=0}^\infty \frac{(-1)^n}{(2n+1)^2}
= \frac{8K}{\pi}
$$
The same result we obtained using contour integrals.
