You can use residues, let $f(z)=\frac{1}{(1+z^2)(1+r^2z^2)(1+r^4z^2)...}$ this has an infite set of singularities at $z=\pm i 1/r^{n},n\in\Bbb N_0$.
We can see that $\large\int_0^\infty\frac{1}{(1+x^2)(1+r^2x^2)(1+r^4x^2)...}dx=\frac{1}{2}\int_{-\infty}^\infty\frac{1}{(1+x^2)(1+r^2x^2)(1+r^4x^2)...}dx$
So we will consider the upper semi circle contour, spanning over $-\infty$ to $+\infty$, only the positive singularities are in this contour, i.e. $z=i\frac{1}{r}$.
Now $\operatorname{Res}(f(z),z=\frac{i}{r^{n}})=\large\lim_{z\to i\frac{1}{r^{n}}}(z-i\frac{1}{r^{n}})f(z)=\large\lim_{z\to i\frac{1}{r^{n}}}\frac{1}{r^n}(r^nz-i)f(z)$
$=\lim_{z\to i\frac{1}{r^{n}}}\frac{1}{r^n(r^nz+i)}\prod_{j=0,j\ne n}^\infty\frac{1}{(1+x^{2}r^{2j})}$
$=\large\frac{1}{2ir^n}\prod_{j=0,j\ne n}^\infty\frac{1}{(1-r^{2j-2n})}$
Now $\frac{1}{2}\int_{-\infty}^\infty\frac{1}{(1+x^2)(1+r^2x^2)(1+r^4x^2)...}dx=\frac{1}{2}2\pi i\sum \operatorname{Res}(f)=\frac{1}{2}\pi i\sum_{n=0}^\infty\frac{1}{ir^n}\prod_{j=0,j\ne n}^\infty\frac{1}{(1-r^{2j-2n})}$
From here I'm not sure how to reach the closed form, but hopefully this helps to show a different approach, even if you are not too well versed in integrals by residue :)
By De Moivre's formula $\sin(x)=\frac{e^{ix}-e^{-ix}}{2i}$ we have the following Fourier sine series:
$$\frac{\sin(3x)\sin(4x)\sin(5x)\cos(6x)}{\sin^2(x)}\\= -\frac{1}{2} \sin(2x)-\frac{1}{2}\sin(4x)+\sin(8x)+\frac{3}{2}\sin(10x)+\frac{3}{2}\sin(12x)+\sin(14x)+\frac{1}{2}\sin(16 x)$$
and:
$$I(n)=\int_{0}^{+\infty}\frac{\sin(2nx)}{x\cosh(x)}\,dx = 2\arctan\left(\tanh\frac{\pi n}{2}\right) $$
follows by differentiation under the integral sign. The original integral can so be expressed in terms of the Gudermannian function:
$$ I = \frac{1}{2} \big(-\text{gd}(\pi)- \text{gd}(2\pi) +
2 \text{gd}(4\pi) + 3 \text{gd}(5\pi) +
3 \text{gd}(6\pi) + 2 \text{gd}(7\pi) +
\text{gd}(8\pi)\big) \approx 7.11363 $$
Best Answer
Proof :
$$\begin{align*} & \displaystyle \sum\limits_{k=0}^{\infty}(-x)^{k}\zeta_H(k+1,a)\\ & \displaystyle = \sum\limits_{k,n=0}^{\infty} \dfrac{(-x)^k}{(n+a)^{k+1}}\\ & \displaystyle = \sum\limits_{n=0}^\infty \dfrac{1}{n+a}\sum\limits_{k=0}^\infty \left(\dfrac{-x}{n+a}\right)^k \\ & \displaystyle = \sum\limits_{n=0}^\infty \dfrac{1}{n+a+x} \\ & \displaystyle = -\psi(a+x)\end{align*}$$
Now integrating,
$$\displaystyle \sum\limits_{k=1}^{\infty} \dfrac{(-1)^{k-1}}{k}\zeta_H(k,a)x^k =\ln\left(\dfrac{\Gamma(a)}{\Gamma(a+x)}\right)\\ $$
Proof :
It is sufficient to evaluate the series,
$$\begin{align*} & \displaystyle\sum\limits_{n=0}^{\infty}\ln\left(1+\dfrac{x}{n+a}\right) \\ & = \displaystyle \sum\limits_{k=1}^\infty\sum\limits_{n=0}^\infty \dfrac{(-1)^{k-1}}{k}x^k \dfrac{1}{(n+a)^k} \\ & =\displaystyle \sum\limits_{k=1}^{\infty} \dfrac{(-1)^{k-1}}{k}\zeta_H(k,a)x^k \\ & =\displaystyle \ln\left(\dfrac{\Gamma(a)}{\Gamma(a+x)}\right)\end{align*}$$
$$\displaystyle \prod\limits_{n=0}^{\infty} \left(1+\dfrac{x}{n+a}\right) = \dfrac{\Gamma(a)}{\Gamma(a+x)}$$
Therefore,
$$\displaystyle \prod\limits_{n=0}^{\infty} \left(1+\dfrac{x^2}{(n+a)^2}\right)\;=\; \dfrac{\Gamma^2 (a)}{\Gamma(a+ix)\Gamma(a-ix)}$$
Proof :
$$\displaystyle F(it)=\int\limits_0^\infty x^{it}\dfrac{f(x)}{x}\; dx$$
Set $x=e^y$ ,
$$\displaystyle F(it)=\int\limits_0^\infty e^{ixt}f(e^x)\; dx$$
Now by properties of Fourier Transform,
$$\begin{align*} & \displaystyle f(e^t)=\int\limits_{-\infty}^\infty g(x)e^{-ixt}\; dx \\ & \displaystyle g(t)=\dfrac{1}{2\pi}\int\limits_{-\infty}^\infty f(e^x)e^{ixt}\; dx\\\end{align*}$$
$$\begin{align*}\displaystyle F(it) & =2\pi g(t) \\ \displaystyle \int\limits_{-\infty}^\infty |F(it)|^2\; dt & = 4\pi^2 \int\limits_{-\infty}^\infty |g(t)|^2\; dt \\ \displaystyle \int\limits_{-\infty}^\infty |F(it)|^2\; dt & = 2\pi \int\limits_{-\infty}^\infty g(t) \int\limits_{-\infty}^\infty e^{ixt}f(e^x)\; dx\; dt \\ \displaystyle \int\limits_{-\infty}^\infty |F(it)|^2\; dt & = 2\pi \int\limits_{-\infty}^\infty f(e^x) \int\limits_{-\infty}^\infty e^{ixt}g(t)\; dt\; dx \\ \displaystyle \int\limits_{-\infty}^\infty |F(it)|^2\; dt & = 2\pi \int\limits_{-\infty}^\infty f(e^x)\overline{f(e^x)}\; dx =2\pi\int\limits_{-\infty}^\infty |f(e^x)|^2\; dx\end{align*}$$
Now by setting $e^x=t$ we get our result,
$$\displaystyle \int\limits_{-\infty}^\infty |F(it)|^2\; dt = 2\pi\int\limits_{-\infty}^\infty \dfrac{|f(t)|^2}{t}\; dt$$
Proof : If we denote the integral by $I$ then using $(2)$ it can be rewritten as,
$$\displaystyle I = \dfrac{\Gamma^2 (b+1)}{\Gamma^2 (a)}\dfrac{1}{2} \int\limits_{-\infty}^\infty \dfrac{|\Gamma(a+ix)|^2}{|\Gamma(b+1+ix)|^2}\; dx$$
Now by defining $\displaystyle h(x)=\dfrac{x^a (1-x)^{b-a}}{\Gamma(b-a+1)}$ for $x\in[0,1]$ and $0$ for $\forall x\notin[0,1]$ we can conclude that $\displaystyle F(s)=M[h(x)]=\dfrac{\Gamma(s+a)}{\Gamma(s+b+1)}$ and from $(3)$ it follows that,
$$\displaystyle I = \dfrac{\Gamma^2 (b+1)}{\Gamma^2 (a)}\dfrac{1}{2} \int\limits_{0}^1 \dfrac{|h(x)|^2}{x}\; dx = \dfrac {\sqrt{\pi}}2\dfrac {\Gamma\left(a+\frac 12\right)\Gamma\left(b+1\right)\Gamma(b-a+1)}{\Gamma(a)\Gamma\left(b+\frac 12\right)\Gamma\left(b-a+\frac 12\right)}$$
where last line follows from the Duplication formula , and we are done !
$$\large \color{red}{\color{blue}{\boxed{\mathfrak{PROVED}}}} $$