First of all, this is not an area:
$$ \int _{-\infty}^\infty e^{-x^{2}}\,dx . $$
It is a definite integral that evaluates to a number. You might associate the number with an area, but only under specific interpretations of what the integral is for. (For example, if someone asks for the area of the region below the graph of the function $e^{-x^{2}}$ and above the $x$-axis, then this integral is how you compute that area.)
In this exercise the fact that there are such interpretations is mostly spurious and misleading.
Just to be clear, you really should not get hung up on having to imagine an area every time you see an integral. There are many applications
(for example, electromagnetism)
that are chock full of integrals, most of which are only related to any kind of "area" in the most abstract way.
So when you see
$$ \int_{-\infty}^\infty e^{-x^{2}}\,dx \int_{-\infty}^\infty e^{-y^{2}}\,dy, $$
it's just the product of two numbers, nothing more.
Certainly not the product of two areas; that makes no sense (at least in this context).
The "two" numbers in fact are just one number:
$$ \int_{-\infty}^\infty e^{-x^{2}}\,dx = \int_{-\infty}^\infty e^{-y^{2}}\,dy.$$
Whether we write $x$ or $y$ or $\theta$ inside the integral really doesn't matter (yet).
It's when we write
$$ \int_{-\infty}^\infty \int_{-\infty}^\infty e^{-(x^2+y^2)}\,dx\,dy $$
that it finally makes sense to regard the $x$ and $y$ as $x$ and $y$ coordinates on a Cartesian plane, because that interpretation helps in visualizing the transformation to polar coordinates (if that's what comes next) or in reinterpreting the volume under $e^{-(x^2+y^2)}$ as a set of concentric shells or stacked disks.
So the progression is number, number times itself (i.e. number squared),
number times itself written slightly differently,
rearrange the order of integration so we have a double integral instead of two single ones, and then
reinterpret the double integral as a volume between a curved surface and the Cartesian plane.
The last step should be the first time we invoke any geometric intuition.
The step
$$ \int_{-\infty}^\infty e^{-x^2}\,dx \int_{-\infty}^\infty e^{-y^2}\,dy =
\int_{-\infty}^\infty \int_{-\infty}^\infty e^{-(x^2+y^2)}\,dx\,dy $$
is admittedly a bit much to swallow all at once if you're not used to this sort of thing. Step by step, it can be done like this:
$$ \int_{-\infty}^\infty e^{-x^2}\,dx = K. $$
This requires that we know (or assume) the integral exists.
But if it exists, it's just a number, so we can call it $K.$
$$ \int_{-\infty}^\infty e^{-x^2}\,dx \int_{-\infty}^\infty e^{-y^2}\,dy
= K \int_{-\infty}^\infty e^{-y^2}\,dy .$$
Simple substitution using the previous equation.
$$ K \int_{-\infty}^\infty e^{-y^2}\,dy = \int_{-\infty}^\infty Ke^{-y^2}\,dy.$$
Since $K$ is simply a fixed number (though we haven't computed it yet), we get the same result multiplying by it before or after integrating.
$$\int_{-\infty}^\infty Ke^{-y^2}\,dy =
\int_{-\infty}^\infty \left(\int_{-\infty}^\infty e^{-x^2}\,dx\right)e^{-y^2}\,dy.
$$
All we did here was reverse the substitution we did earlier.
Now we look at what is inside the $dy$ integral:
$$ \left(\int_{-\infty}^\infty e^{-x^2}\,dx\right)e^{-y^2}.$$
In this expression $e^{-y^2}$ is simply an (unknown) number that we are multiplying the integeral $\int_{-\infty}^\infty e^{-x^2}\,dx$ by, so again we can get the same result multiplying inside the integral as outside:
$$ e^{-y^2}\int_{-\infty}^\infty e^{-x^2}\,dx =
\int_{-\infty}^\infty e^{-y^2}e^{-x^2}\,dx.$$
Finally, $e^{-y^2}e^{-x^2} = e^{-(x^2+y^2)}.$
Putting it all together we get
$$
\int_{-\infty}^\infty \left(\int_{-\infty}^\infty e^{-x^2}\,dx\right)e^{-y^2}\,dy
= \int_{-\infty}^\infty \left(\int_{-\infty}^\infty e^{-(x^2+y^2)}\,dx\right)dy.
$$
And then Fubini's Theorem lets us treat
$$
\int_{-\infty}^\infty \left(\int_{-\infty}^\infty e^{-(x^2+y^2)}\,dx\right)dy
= \int_{-\infty}^\infty \int_{-\infty}^\infty e^{-(x^2+y^2)}\,dx\,dy
$$
as an integral over the combined $(x,y)$ coordinates of a Cartesian plane
instead of just one integral inside another,
which we need in order to find the circles in the shape of this integral.
Conceptually, the direction in which I would approach the calculation of the constant in the Gaussian distribution is exactly the opposite of the direction taken in the previous part of this answer.
That is, I would start by setting up a joint distribution of two iid Gaussian variables as a function over an $x,y$ Cartesian plane,
using an as-yet-unknown constant factor to make this a probability distribution.
That is, I would start with a distribution that is already two-dimensional.
Then I would show that the integral of that distribution can be written as the product of two integrals, which can be rewritten as the square of one integral.