Here is an elementary proof. Using the formula
\begin{equation}
\frac1{a^n}=\frac1{\Gamma(n)}\,\int_0^\infty u^{n-1} e^{-u\,a}\,du \tag{1}
\end{equation}
with $a=1+\|x\|^2$,
denoting your integral by $J_n$, letting $I:=[0,1]$, $\phi(z):=\frac1{\sqrt{2\pi}}\,e^{-z^2/2}$, and $\Phi(z):=\int_{-\infty}^z\phi(t)\,dt$, and making substitutions $x_1=z_1/\sqrt{2u}$ and then $u=z^2/2$, we have
\begin{align}
J_n&=\frac1{\Gamma(n)}\,\int_0^\infty u^{n-1} e^{-u}\,du \int_{I^{2n-1}}e^{-u\,\|x\|^2}\,dx \\
&=\frac1{\Gamma(n)}\,\int_0^\infty u^{n-1} e^{-u}\,du \Big(\int_I e^{-u\,x_1^2}\,dx_1\Big)^{2n-1} \\
&=\frac{\pi^{n-1/2}}{\Gamma(n)}\,\int_0^\infty u^{-1/2} e^{-u}\,du \Big(\Phi(\sqrt{2u})-\frac12\Big)^{2n-1} \\
&=\frac{2\pi^{n}}{\Gamma(n)}\,\int_0^\infty \Big(\Phi(z)-\frac12\Big)^{2n-1}\,\phi(z)\,dz \\
&=\frac{2\pi^{n}}{\Gamma(n)}\,\int_0^\infty \Big(\Phi(z)-\frac12\Big)^{2n-1}\,d\Phi(z) \\
&=\frac{\pi^{n}}{4^n \Gamma(n+1)},
\end{align}
as desired.

This derivation obviously holds whenever $2n-1$ is a positive integer.

**Further comments:**

More generally, in view of Bernstein's theorem on completely monotone functions, one can quite similarly express the box integral
\begin{equation}
\int_{I^n} g(\|x\|^2)\,dx
\end{equation}
as an ordinary integral over $[0,\infty)$ -- for any completely monotone function $g$. Indeed, Bernstein's theorem states that any completely monotone function is a positive mixture of decreasing exponential functions; identity (1) is then such an explicit Bernstein representation of the completely monotone function $a\mapsto \frac1{a^n}$.

Now, furthermore, one does not have to insist that the mixture be positive, since we are dealing here with identities rather than inequalities. So, any function $g$ that can be represented as the mixture
\begin{equation}
g(a)=\int_0^\infty e^{-u\,a}\,\mu(du)
\end{equation}
of decreasing exponential functions $a\mapsto e^{-u\,a}$
for $a>0$, where $\mu$ is any finite, possibly signed ("mixing") measure, will do just as well.