In this case you don't need Fubini's theorem or equivalent to produce an iterated integral (although a careless exposition can conceal this fact), since what one is actually doing is using the linearity of the integral twice: you have
$$ J^2 = J\int_0^{\infty} e^{-x^2} \, dx = \int_0^{\infty} e^{-x^2} J \, dx. $$
Now, $e^{-x^2}$ is also a constant when integrating with respect to $y$, so we have
$$ e^{-x^2} J = e^{-x^2} \int_0^{\infty} e^{-y^2}\, dy = \int_0^{\infty} e^{-x^2}e^{-y^2} \, dy, $$
and so
$$ J^2 = \int_0^{\infty} \left( \int_0^{\infty} e^{-(x^2+y^2)} dy \right) dx $$
as desired.
The inner integral $\int_0^{\infty} e^{-(x^2+y^2)} dy$ is now a function of $x$, and we substitute to write this function in a different way, as $x\int_0^{\infty} e^{-x^2(1+t^2)} \, dt $. Now is when you are forced to change the order of integration using some theorem or other. Either we embrace a Lebesgue integral version, such as Tonelli's theorem,
Let $f : (a,b) \times (c,d) \to \mathbb{R}$ be positive and Lebesgue-integrable, where $-\infty \leq a,b \leq \infty$ and $ -\infty \leq c,d \leq \infty$. Then
$$ \int_{(a,b)\times (c,d)} f(x,y) \, dx \times dy = \int_a^b \left( \int_c^d f(x,y) \, dy \right) dx = \int_c^d \left( \int_a^b f(x,y) \, dx \right) dy. $$
Here, though, we can manage with a simple version for Riemann integrals:
Let $f : (a,b) \times (c,d) \to \mathbb{R}$ be positive and Riemann-integrable, where $-\infty <a,b < \infty$ and $ -\infty < c,d < \infty$. Then
$$ \int_{(a,b)\times (c,d)} f(x,y) \, dx \times dy = \int_a^b \left( \int_c^d f(x,y) \, dy \right) dx = \int_c^d \left( \int_a^b f(x,y) \, dx \right) dy. $$
In particular, $x = 1-(1-u)^{-1}$ provides a bijection between $(0,1)$ and $(0,\infty)$, so we find that
$$ J^2 = \int_0^1 \int_0^1 \frac{u}{(1-u)^3(1-v)^2}\exp{\left( - \frac{u^2}{(1-u)^2}\left( 1 + \frac{v^2}{(1-v)^2} \right) \right)} \, dv \, du. $$
It turns out that the integrand is in fact continuous on the whole open square, so it's integrable, and we can use the lightweight Tonelli to swap the integrals, and then after undoing the substitutions, we're finally at
$$ J^2 = \int_0^{\infty} \int_0^{\infty} xe^{-x^2(1+t^2)} \, dx \, dt, $$
which is straightforward to do. Alternatively if one is happy to use the Lebesgue integral, we don't need to muck about with substitutions, and can apply the proper version of Tonelli to the integral on $(0,\infty) \times (0,\infty)$.
Suppose $\mu_1=\mu_2$ are counting measures on $\Omega_1=\Omega_2=\{1,2,\ldots\}$.
Define the following function on $\Omega_1\times\Omega_2$:
$$f(i,j)=\begin{cases}1&,\text{ if }i=j
\\ -1&,\text{ if }i=j+1
\\ 0&,\text{ otherwise }
\end{cases}$$
We can write out the values of $f(i,j)$ in a matrix form like
$$[f(i,j)]=\begin{bmatrix}1&0&0&0&\cdots
\\ -1&1&0&0&\cdots
\\0&-1&1&0&\cdots
\\0&0&-1&1&\cdots
\\\vdots&\vdots&0&-1&\cdots
\\\vdots&\vdots&\vdots&\vdots&\ddots
\\0&0&0&0&\cdots
\end{bmatrix}$$
Only the first row sums to $1$, each of the remaining rows sum to $0$. Also sum of each column is $0$.
Therefore, $$\int\left(\int f(x,y)\,d\mu_2(y)\right)d\mu_1(x)=\sum_{i=1}^\infty \left(\sum_{j=1}^\infty f(i,j)\right)=1$$
And $$\int\left(\int f(x,y)\,d\mu_1(x)\right)d\mu_2(y)=\sum_{j=1}^\infty \left(\sum_{i=1}^\infty f(i,j)\right)=0$$
However,
\begin{align}
\iint|f(x,y)|\,d\mu_1(x)\,d\mu_2(y)&=\sum_{i=1}^\infty\sum_{j=1}^\infty|f(i,j)|
\\&=\sum_{i=1}^\infty\left(\sum_{j=1}^\infty |f(i,j)|\right)\quad,\small\text{ by Fubini/Tonelli, since }|f|\ge 0
\\&=1+2+2+\cdots
\\&=\infty
\end{align}
So $f$ is not $\mu$-integrable where $\mu=\mu_1\otimes\mu_2$ is the product measure.
Best Answer
If we have Tonelli's theorem, then by considering positive parts and negative parts separately we immediately obtain Fubini's theorem.
Conversely, assuming Fubini's theorem, Tonelli's theorem follows by monotone convergence argument applied to cut-off functions $f_k(x) = \min \{k, f(x)\} \chi_{B_k}(x)$. You can also find the detail at the Chapter 6.2 of the celebrated textbook Measure and Integral by Wheeden and Zygmund.