Here is a way you can do it with brute force with the aid of change of variables and Fubini's Theorem. Let
\begin{align*}
I &= \int_{0}^\infty \frac{1}{(1+x^2)(1+r^2x^2) \ \dots (1+r^{2n}x^2)} \ dx.
\end{align*}
Note for each $i \in \lbrace 1, \ \dots \ ,n \rbrace,$ we can write
\begin{align*}
\frac{1}{1+r^{2i} x^2} &= \int_{0}^{\infty} e^{-(1+r^{2i}x^2) y_i} \ dy_i.
\end{align*}
Thus, using Fubini's Theorem,
\begin{align*}
I &=\int_{0}^\infty \int_{0}^\infty \ \dots \ \int_{0}^\infty e^{-y_1(1+x^2)} e^{-y_2(1+r^2x^2)} \ \dots \ e^{-y_{n}(1+r^{2n}x^2)} \ dy_n \ \dots \ dy_1\ dx \\
&=\int_{0}^\infty \int_{0}^\infty \ \dots \ \int_{0}^\infty e^{-(y_1 + \ \dots \ + y_n)} \ e^{-(y_1+r^2y_2 \ \dots \ + r^{2n}y_n)x^2} \ dy_n \ \dots \ dy_1 \ dx \\
&=\int_{0}^\infty \int_{0}^\infty \ \dots \ \int_{0}^\infty e^{-(y_1 + \ \dots \ + y_n)} \ e^{-(y_1+r^2y_2 \ \dots \ + r^{2n}y_n)x^2} \ dx \ dy_n \ \dots \ dy_1.
\end{align*}
To carry out the integral with respect to $x,$ use the one dimensional change of variables $x=\frac{t}{\sqrt{y_1+r^2y_2 \ \dots \ + r^{2n}y_n}}$ and the well-known Gaussian integral $\int_{0}^\infty e^{-t^2} \ dt= \frac{\sqrt{\pi}}{2}.$ We get
\begin{align*}
I &= \frac{\sqrt{\pi}}{2} \int_{0}^{\infty} \ \dots \ \int_{0}^{\infty} \frac{e^{-(y_1 + \ \dots \ + y_n)}}{\sqrt{y_1+r^2y_2 \ \dots \ + r^{2n}y_n}} \ dy_n \ \dots \ dy_1.
\end{align*}
Now, perform the multidimensional change of variables
\begin{align*}
y_i &= t_i - r^2 t_{i+1}, \quad 1 \leq i <n \\
y_n &= t_n.
\end{align*}
By induction on $n,$ it is easy to verify $\left|\frac{\partial(y_1 ,\ \dots \ , y_n)}{\partial(t_1 , \ \dots \ , t_n)}\right|=1$ and that transformed region of integration is
\begin{align*}
t_1 &>0 \\
0 &< t_i < \frac{t_{i-1}}{r^2},\quad 2 \leq i \leq n.
\end{align*}
With this,
\begin{align*}
I &= \frac{\sqrt{\pi}}{2} \int_{0}^{\infty} \int_{0}^{t_1/r^2} \ \dots \ \int_{0}^{t_{n-1}/r^2} \frac{e^{-t_1-(1-r^2)(t_2 + \ \dots \ + t_n)}}{\sqrt{t_1}}\ dt_n \ \dots \ dt_2 \ dt_1.
\end{align*}
This makes $I$ an "elementary" integral that can be evaluated either by hand or with Mathematica. To carry out the integration with respect to $t_n, \ \dots \ , t_2,$ one must repeatedly apply the exponential identity
$$\int_{0}^{t} e^{-ax} \ dx = \frac{1-e^{-at}}{a}$$ for $a>0,$ which is pretty messy. To carry out the final integration with respect to $t_1,$ make the substitution $t_1=u_1^2$ which will transform the resulting integrand into a Gaussian function.
At the moment, I do not know of a clean, slick way to evaluate this multidimensional integral and arrive at the general answer, but I will look into it.