For ease of exposition, let us consider zero-mean unit-variance jointly continuous random variables $X$ and $Y$ with joint density $f(x,y)$. Then,
$$\begin{align}
\rho = E[XY] &= \int_{-\infty}^\infty \int_{-\infty}^\infty
x\cdot y\cdot f(x,y)\,\mathrm dx\,\mathrm dy\\
&= \int_0^\infty \int_0^\infty xy\big[f(x,y)+f(-x-y)-f(-x,y)-f(x,-y)\big]
\,\mathrm dx\,\mathrm dy.
\end{align}$$
If $\rho > 0$, then the term
$\big[f(x,y)+f(-x-y)-f(-x,y)-f(x,-y)\big]$ cannot be negative for all
$(x,y)$ in the first quadrant; indeed, roughly speaking,
$$f(x,y)+f(-x-y) > f(-x,y)+f(x,-y)$$
quite often, that is, there is generally more mass in the vicinity of $(x,y)$ and $(-x-y)$ (i.e. in the first and third quadrants) than in the
vicinity of $(-x,y)$ and $(x,-y)$, (i.e. in the second and fourth quadrants). For $\rho$ close to $1$, very roughly speaking, the probability mass lies in the vicinity of the line $x=y$.
For $\rho < 0$ and for $\rho$ close to $-1$, similar remarks apply
mutatis mutandis.
The simulations in
this answer of Michael Hardy illustrate this notion very well. My point is that this sort
of behavior (more mass in two opposite quadrants than in the complementary
opposite quadrants when $|\rho|$ is close to $1$) is a general feature shared by all (zero-mean unit-variance) random variables. Similar
remarks apply when the mean point is not the origin (the quadrants
are defined w.r.t. the mean point) and the variances are different
(the slopes of the lines near which the mass lies are different)