The sum of two independent brownian motions is also a brownian motion (which is formally proved here)
First of all we solve for $Y_t$. We notice that $d\bigl( Y e^{\lambda t} \bigr) = dY e^{\lambda t} + \lambda Y e^{\lambda t}$.
\begin{align}
dY_t &= -\lambda Y_t dt + \sigma_Y \rho dZ_{1t} + \sigma_Y \sqrt{1 - \rho^2} dZ_{2t} \\
&= -\lambda Y_t dt + \sigma_Y dW_t \\
dY_t + \lambda Y dt &= \sigma_Y dW_t \\
d\bigl( Y_t e^{\lambda t} \bigr) &= e^{\lambda t} \sigma_Y dW_t \\
\Bigl[ Y_t e^{\lambda t} \Bigr]^T_0 &= \sigma_Y \int^T_0 e^{\lambda t} dW_t \\
Y_T &= \sigma_Y e^{-\lambda T}\int^T_0 e^{\lambda t} dW_t \\
\end{align}
Where I've defined $W_t = \rho Z_{1t} + \sqrt{1 - \rho^2} Z_{2t}$. This is gaussian distributed with an expectation of 0 and a variance ${\frac {\sigma^2_Y} {2 \lambda}} \bigl( 1 - e^{-2\lambda T} \bigr)$ coming from the Ito isometry
Now plugging this in and solving for $X_t$:
\begin{align}
dX_t &= Y_t dt + \sigma_X dZ_{1t} \\
&= \bigl( \sigma_Y e^{-\lambda t}\int^t_0 e^{\lambda s} dW_s \bigr) dt + \sigma_X dZ_{1t} \\
\Bigl[ X_t \Bigr]^T_0 &= \sigma_Y \int^T_0 e^{-\lambda t} \bigl( \int^t_0 e^{\lambda s} dW_s \bigr) dt + \sigma_X Z_{1t}
\end{align}
We can solve $\int^T_0 e^{-\lambda t} \bigl( \int^t_0 e^{\lambda s} dW_s \bigr) dt$ using stochastic integration by parts, as done here, using $A = \int^t_0 e^{-\lambda s} ds$ and $B = \int^t_0 e^{\lambda s} dW_s$ gives
\begin{align}
\Bigl[ A_t \cdot B_t \Bigr]^T_0 &= \int^T_0 e^{-\lambda t} \bigl( \int^t_0 e^{\lambda s} dW_s \bigr) dt + \int^T_0 \bigr( \int^t_0 e^{-\lambda s} ds \bigr) e^{\lambda t} dW_t \\
\int^T_0 e^{-\lambda t} \bigl( \int^t_0 e^{\lambda s} dW_s \bigr) dt &= -\int^T_0 \bigr( \int^t_0 e^{-\lambda s} ds \bigr) e^{\lambda t} dW_t + \bigl( \int^T_0 e^{-\lambda s} ds \bigr) \cdot \bigl( \int^T_0 e^{\lambda t} dW_t \bigr) \\
&= -{\frac {1} \lambda}\int^T_0 (e^{\lambda t} - 1) dW_t + {\frac {1} \lambda} (1 - e^{-\lambda T}) \int^T_0 e^{\lambda t} dW_t \\
&= -{\frac {1} \lambda}\int^T_0 (e^{\lambda t} - 1) - e^{\lambda t}(1 - e^{-\lambda T}) dW_t\\
&= {\frac {1} \lambda}\int^T_0 \bigl(1 - e^{-\lambda (T-t)} \bigr) dW_t
\end{align}
And substituting this in above, we have
\begin{align}
X_T &= {\frac {\sigma_Y} \lambda}\int^T_0 \bigl(1 - e^{-\lambda (T-t)} \bigr) dW_t + \sigma_X Z_{1T}
\end{align}
This is the sum of two (correlated) gaussians, so it is also a gaussian as required
- From above, we have
\begin{align}
\Bigl[ Y_s e^{\lambda s} \Bigr]^{T+t}_t &= \sigma_Y \int^{T+t}_t e^{\lambda s} dW_s \\
Y_{T+t} &= e^{-\lambda T} Y_t + e^{-\lambda (T+t)} \sigma_Y \int^{T+t}_t e^{\lambda s} dW_s
\end{align}
Conditioning on $Y_t$, we can now find $X_{T+t}$ as above
\begin{align}
{\mathbb E}\bigl[(X_{T+t} - X_t \bigr)] &= {\mathbb E}\bigl[ \int_t^{T+t} dX_s \bigr] \\
&= {\mathbb E}\bigl[\int^{T+t}_t Y_s ds + \int^{T+t}_t \sigma_X dZ_{1s} \bigr]\\
&= {\mathbb E}\bigl[\int^{T+t}_t Y_s ds\bigr]\\
&= {\mathbb E}\bigl[\int^{T}_0 Y_{u+t} du\bigr] \\
&= {\mathbb E}\bigl[\int^{T}_0 \Bigl( e^{-\lambda u} Y_t + e^{-\lambda (u+t)} \sigma_Y \int^{u+t}_t e^{\lambda s} dW_s \Bigr) du \bigr] \\
&= {\mathbb E}\bigl[\int^{T}_0 e^{-\lambda u} Y_t du \bigr]\\
&= {\frac 1 {\lambda}} Y_t \bigl( 1 - e^{-\lambda T} \bigr)
\end{align}
(where I've changed variables from $s$ to $u = s - t$) which makes sense - $Y_t$ is mean-reverting so we expect future values to be closer to zero than current values
\begin{align}
{\mathbb E}\bigl[(X_{T+t} - X_t \bigr)^2] &= {\mathbb E}\bigl[ \bigl( \int_t^{T+t} dX_s \bigr)^2 \bigr] \\
&= {\mathbb E}\bigl[\Bigl(\int^{T+t}_t Y_s ds + \int^{T+t}_t \sigma_X dZ_{s1} \Bigr)^2 \bigr]\\
&= {\mathbb E}\bigl[\bigl( \int^{T+t}_t Y_s ds \bigr)^2 + \int^{T+t}_t \sigma_X^2 dt + 2 \int^{T+t}_t Y_s ds \int^{T+t}_t \sigma_X dZ_{1s} \bigr]\\
&= {\frac {Y_t^2} {\lambda^2}} \bigl( 1 - e^{-\lambda T} \bigr)^2 + \sigma_X^2 T + 2 {\mathbb E}\bigl[ \int^{T+t}_t Y_s ds \int^{T+t}_t \sigma_X dZ_{1s} \bigr]
\end{align}
For clarity I break out the last term separately:
\begin{align}
{\mathbb E}\bigl[ \int^{T+t}_t Y_s ds \int^{T+t}_t \sigma_X dZ_{1s} \bigr] &= {\mathbb E}\bigl[ \int^T_0 Y_{u+t} du \int^T_0 \sigma_X dZ_{u1} \bigr] \\
&= {\mathbb E}\bigl[ \int^T_0 \Bigl( e^{-\lambda u} Y_t + e^{-\lambda (u+t)} \sigma_Y \int^{u+t}_t e^{\lambda s} dW_s \Bigr) du \int^T_0 \sigma_X dZ_{1u} \bigr] \\
&= {\mathbb E}\bigl[ \int^T_0 \Bigl( e^{-\lambda (u+t)} \sigma_Y \int^{u+t}_t e^{\lambda s} dW_s \Bigr) du \int^T_0 \sigma_X dZ_{1u} \bigr] \\
&= {\mathbb E}\bigl[ \int^T_0 \Bigl( e^{-\lambda (u+t)} \sigma_Y \int^{u+t}_t e^{\lambda s} \rho dZ_{1u} \Bigr) du \int^T_0 \sigma_X dZ_{1u} \bigr] \\
&= {\frac {\rho \sigma_Y} {\lambda}} {\mathbb E}\bigl[ \int^T_0 \bigl( 1 - e^{-\lambda T} \bigr) dZ_{1u} \int^T_0 \sigma_X dZ_{1u} \bigr] \\
&= {\frac {\rho \sigma_X \sigma_Y} {\lambda}} \int^T_0 \bigl( 1 - e^{-\lambda T} \bigr) du\\
&= {\frac {T \rho \sigma_X \sigma_Y} {\lambda}}\bigl( 1 - e^{-\lambda T} \bigr)
\end{align}
and plugging this back in to the block above we have
\begin{align}
{\mathbb E}\bigl[(X_{T+t} - X_t \bigr)^2] &= {\frac {Y_t^2} {\lambda^2}} \bigl( 1 - e^{-\lambda T} \bigr)^2 + \sigma_X^2 T + 2 {\frac {T \rho \sigma_X \sigma_Y} {\lambda}} \bigl( 1 - e^{-\lambda T} \bigr)
\end{align}
Thinking about the behaviour of this process as $T \to \infty$, we see that the $( 1 - e^{-\lambda T})$ terms go to zero and we're left with a variance of $\sigma_X^2 T$, which is just standard geometric brownian motion's variance.
As $T \to 0$, $( 1 - e^{-\lambda T}) \to \lambda T$ which cancel out all of the $\lambda$s so the expression becomes
\begin{align}
{\frac 1 T} \lim_{T \to 0} {\mathbb E}\bigl[(X_{T+t} - X_t \bigr)^2] &= Y_t^2 T + \sigma_X^2 - 2 \rho \sigma_X \sigma_Y T
\end{align}
so as expected, short term variance at $t$ increases with the level of $Y_t$, and decreases if the processes are more positively correlated.
Wow what a question!