Rather than embarking on some pretty involved computations of conditional distributions, one should rely on one of the main assets of Gaussian families, namely, the...
Key feature: In Gaussian families, conditioning acts as a linear projection.
Hence, as the OP suggested, one could do worse than to start from a representation of $(X,Y)$ by standard i.i.d. Gaussian random variables $U$ and $V$, for example,
$$
X=\mu_x+\sigma_xU\qquad
Y=\mu_y+\sigma_y(\rho U+\tau V)$$ where the parameter $\tau$ is $$\tau=\sqrt{1-\rho^2}
$$
Since $\sigma_x\ne0$, the sigma-algebra generated by $X$ is also the sigma-algebra generated by $U$ hence conditioning by $X$ or by $U$ is the same. Furthermore, constants and functions of $X$ or $U$ are all $U$-measurable while functions of $V$ are independent on $U$, thus,
$$
\mathrm E(Y\mid X)=\mu_y+\sigma_y(\rho U+\tau \mathrm E(V))=\mu_y+\sigma_y \rho U
$$
which is equivalent to
$$
\color{red}{\mathrm E(Y\mid X)=\mu_y+\rho\frac{\sigma_y}{\sigma_x}(X-\mu_x)}
$$
Likewise, when computing conditional variances conditionally on $X$, deterministic functions of $X$ or $U$ should be considered as constants, hence their conditional variance is zero, and functions of $V$ are independent on $X$, hence their conditional variance is their variance. Thus,
$$
\mbox{Var}(Y\mid X)=\mbox{Var}(\sigma_y\tau V\mid X)=\sigma_y^2\tau^2\mbox{Var}(V\mid X)=\sigma_y^2\tau^2\mbox{Var}(V)
$$
that is,
$$
\color{red}{\mbox{Var}(Y\mid X)=\sigma_y^2(1-\rho^2)}
$$
Finally, the event $$A=[X>\mu_x,Y>\mu_y]$$ is also
$$
A=[U>0,\rho U+\tau V>0].
$$
To evaluate $\mathrm P(A)$, one can turn to the planar representation of couples of independent standard Gaussian random variables, which says in particular that the distribution of $(U,V)$ is invariant by rotations. The event $A$ means that the direction of the vector $(U,V)$ is between the angle $\vartheta$ in $(-\pi/2,\pi/2)$ such that $$\tan(\vartheta)=-\rho/\tau$$ and the angle $\pi/2$. Thus,
$$\mathrm P(A)=\frac{\pi/2-\vartheta}{2\pi}$$
that is,
$$
\color{red}{\mathrm P(X>\mu_x,Y>\mu_y)=\frac14+\frac1{2\pi}\arcsin\rho}
$$
Numerical application: If $\mu_x=2$, $\mu_y=-1$, $\sigma_x=2$, $\sigma_y=1$ and $\rho=-\sqrt3/2$, then
$$
\mathrm E(Y\mid X)=-1+\sqrt3/2-(\sqrt3/4)X\qquad
\mbox{Var}(Y\mid X)=1/4
$$
and $\tau=1/2$, hence $\vartheta=\pi/3$ and $$\mathrm P(A)=1/12$$
We have that
$$\begin{align}
\operatorname E\mathrm{sgn}(X-Y)
& =-1\cdot\Pr(X<Y)+0\cdot\Pr(X<Y)+1\cdot\Pr(X>Y)
\\ & =\Pr(X>Y)-\Pr(X<Y)
\end{align}$$
and
$$\begin{align}
\operatorname E\mathrm{sgn}^2(X-Y)
& =(-1)^2\cdot\Pr(X<Y)+0^2\cdot\Pr(X<Y)+1^2\cdot\Pr(X>Y)
\\ & =\Pr(X>Y)+\Pr(X<Y)
\end{align}$$
using the law of the unconscious statistician.
Hence,
$\begin{align}
\operatorname{Var}\mathrm{sgn}(X-Y)
&=\Pr(X>Y)+\Pr(X<Y)-[\Pr(X>Y)-\Pr(X<Y)]^2.
\end{align}$
Best Answer
You really should look at $-Y$ as Normal with mean $-M_y$, variance $\Sigma_y$ and use the formula for $Var (X+Y)$.
Alternatively, for any random variables
$E(X-Y) = EX - EY$ and $Var(X-Y) = Var (X) + Var (Y) - 2 Cov (X,Y)$.
Your expression for Variance is almost always wrong. The only obvious casee I can see it is correct is when $Y$ is a constant or $Y=-X$.
What is interesting here is whether $Z$ is normally distributed. The answer is
Yes if (X,Y) is jointly normal
Not necessarily if (X,Y) are only known to be marginally normal
Consider this construction:
take $X$ to be $N(0,1)$ distribution. Take $Y=X$ if $|X|\leq c$, $Y=-X$, if $|X|>c$, where $c>0$.
Convince yourself $Y$ is also distributed as $N(0,1)$!
here $X-Y=0$ if $|X|\leq c$, else it equals to $2X$, which is never not normal.