$I$ is a number. You obtain $I$ by integrating the Gaussian over the real line. That is,
$$ I := \int_{\mathbb{R}} \mathrm{e}^{-x^2}\,\mathrm{d}x. $$
Notice that after the integration, the variable $x$ is gone. It plays no role, and is just a name for a quantity which varies over $\mathbb{R}$ in the integration. Before going any farther, perhaps it would be useful to consider the following exercise:
Exercise: Compute
$$\int_{0}^{1} x^2\,\mathrm{d}x
\qquad\text{and}\qquad
\int_{0}^{1} y^2\, \mathrm{d}y.
$$
What do you notice about these two quantities?
What you should notice is that the name of the variable with respect to which we are integrating doesn't matter. In this sense, both $x$ and $y$ are "dummy" variables in the integration. They are not elements of $\mathbb{R}^2$ and are not orthogonal to each other. They are simply labels.
So, we know how $I$ is defined, but the name of the variable doesn't matter. Hence we can write
$$ I = \int_{\mathbb{R}} \mathrm{e}^{-y^2}\,\mathrm{d}y $$
just as easily as the original definition. Or even
$$
I
=\int_{\mathbb{R}} \mathrm{e}^{-s^2}\,\mathrm{d}s
=\int_{\mathbb{R}} \mathrm{e}^{-t^2}\,\mathrm{d}t.\tag{1}
$$
The actual names are completely irrelevant. Taking this last version, we can write
\begin{align*}
I^2
&= \left( \int_{\mathbb{R}} \mathrm{e}^{-s^2}\,\mathrm{d}s \right)^2 \\
&= \left( \int_{\mathbb{R}} \mathrm{e}^{-s^2}\,\mathrm{d}s \right)\left( \int_{\mathbb{R}} \mathrm{e}^{-s^2}\,\mathrm{d}s \right) \\
&= \left( \int_{\mathbb{R}} \mathrm{e}^{-s^2}\,\mathrm{d}s \right)\left( \int_{\mathbb{R}} \mathrm{e}^{-t^2}\,\mathrm{d}t \right) && \tag{using (1)} \\
&= \int_{\mathrm{R}} \mathrm{e}^{-s^2} \left( \int_{\mathbb{R}} \mathrm{e}^{-t^2}\,\mathrm{d}t \right) \, \mathrm{d}s \tag{2} \\
&= \int_{\mathrm{R}} \left( \int_{\mathbb{R}} \mathrm{e}^{-s^2} \mathrm{e}^{-t^2}\,\mathrm{d}t\right) \, \mathrm{d}s \tag{3} \\
&= \iint_{\mathbb{R}^2} \mathrm{e}^{-(s^2+t^2)}\,\mathrm{d}(x,y).
\end{align*}
At (2), we are using the fact that
$$ \left( \int \mathrm{e}^{-s^2} \,\mathrm{d}s\right) k
= \int_{\mathbb{R}} \mathrm{e}^{-s^2}k \,\mathrm{d}s $$
for any constant $k$ which does not depend on $s$, and the fact that
$$ \int_{\mathbb{R}} \mathrm{e}^{-t^2} \,\mathrm{d}t $$
is a constant (since it doesn't depend on $s$). At (3), we are playing the same game, but this time using the fact that the value of $\mathrm{e}^{-s^2}$ does not depend on $t$. Once we get the integral at (3), we can apply the Fubini-Tonelli in order to interpret the iterated integral as an integral over a region, then make a change of variables and work with polar (rather than Cartesian) coordinates
To give a simpler example, observe that (via similar arguments)
$$ 1 = \int_{0}^{1} 2s \,\mathrm{d}s. $$
Now, observe that
\begin{align}
1 = 1^2
&= \left( \int_{0}^{1} 2s \,\mathrm{d}s \right)^2 \\
&= \left( \int_{0}^{1} 2s \,\mathrm{d}s \right) \left( \int_{0}^{1} 2t \,\mathrm{d}t \right) \\
&= \int_{0}^{1} 2s \left( \int_{0}^{1} 2t \,\mathrm{d}t \right) \,\mathrm{ds} \\
&= \int_{0}^{1} \left( \int_{0}^{1} (2s)(2t)\,\mathrm{d}t\right)\, \mathrm{d}s \\
&= 4 \int_{0}^{1} \left( \int_{0}^{1} st \,\mathrm{d}t\right)\, \mathrm{d}s \tag{4} \\
&= 4 \int_{0}^{1} \left. \frac{1}{2}st^2 \right|_{t=0}^{1} \,\mathrm{ds} \\
&= 2 \int_{0}^{1} s\, \mathrm{d}s \\
&= 2 \cdot \left.\frac{1}{2}s^2 \right|_{s=0}^{1} \\
&= 1,
\end{align}
which is the expected result. Here, we can evaluate the iterated integral obtained at (4) directly (i.e. no change of variables is required), but otherwise this example follows the same principles as the argument for integrating the Gaussian.