Helping with the essential steps:
Writing
$$\tag{1}
Y_t=\alpha\,\theta+\alpha\sqrt{2\mu}\left(\int_0^t e^{\mu(s-t) }\,dW_s \right)
$$
this SDE is more an inhomogeneous ODE of the form
$$\tag{2}
\dot X_t=-\alpha X_t+Y_t
$$
which has the solution
$$\tag{3}
X_t=\textstyle X_0\,e^{-\alpha t}+\int_0^tY_u\,e^{\alpha(u-t)}\,du\,.
$$
The mean of $Y_t$ is $\mathbb E[Y_t]=\alpha\,\theta\,.$ Therefore, the mean of
$X_t$ is
$$\tag{4}
\mathbb E[X_t]=\textstyle X_0\,e^{-\alpha t}+\int_0^t\alpha\,\theta\,e^{\alpha(u-t)}\,du=X_0\,e^{-\alpha t}+\theta\,(1-e^{-\alpha t})\,.
$$
Further
\begin{align}\tag{5}
\mathbb E[X_t^2]&=\textstyle X_0^2\,e^{-2\alpha t}+2\,X_0\,e^{-\alpha t}\,\mathbb E\Big[\int_0^tY_u\,e^{\alpha(u-t)}\,du\Big]+\mathbb E\Big[\Big(\int_0^tY_u\,e^{\alpha(u-t)}\,du\Big)^2\Big]\,.
\end{align}
The second term in (5) is (similar to the calculation for the mean)
$$\tag{6}
2\,X_0\,e^{-\alpha t}\,\theta\,(1-e^{-\alpha t})=2\,X_0\,\theta\,(e^{-\alpha t}-1)\,.
$$
The last term in (5) is
\begin{align}\tag{7}
\textstyle\mathbb E\Big[\int_0^t\int_0^tY_u\,Y_v\,e^{\alpha(u+v-2t)}\,du\,dv\Big]
=\int_0^t\int_0^t\mathbb E[Y_u\,Y_v]\,e^{\alpha(u+v-2t)}\,du\,dv\,.
\end{align}
Using (1) we get
\begin{align}
\mathbb E[Y_uY_v]&=\alpha^2\,\theta^2\,+\sigma^2\,\mathbb E\Big[\Big(\textstyle\int_0^u e^{\mu(s-u)}\,dW_s\Big)\Big(\textstyle\int_0^v e^{\mu(s-v)}\,dW_s\Big)\Big]\\[2mm]
&=\alpha^2\,\theta^2\,+\sigma^2\,e^{-\mu(u+v)}\mathbb E\Big[\Big(\textstyle\int_0^u e^{\mu s}\,dW_s\Big)\Big(\textstyle\int_0^v e^{\mu s}\,dW_s\Big)\Big]\,.\tag{8}
\end{align}
Because $W_t$ has independent increments and the integrands in the
stochastic integrals are deterministic the last expectation in (8) is
$$\tag{9}
\textstyle\int_0^{\min(u,v)} e^{2\mu s}\,ds=\displaystyle\frac{e^{2\mu\min(u,v)}-1}{2\mu}\,.
$$
(This is a relatively simple generalization of the covariance function $\mathbb E[W_uW_v]=\min(u,v)$ and of the Ito isometry $\mathbb E[(\int_0^ue^{\mu s}\,dW_s)^2]=\int_0^u e^{2\mu s}\,ds\,.$)
Therefore
\begin{align}
\mathbb E[Y_uY_v]
&=\alpha^2\,\theta^2\,+\,\sigma^2\,e^{-\mu(u+v)}\,\frac{e^{2\mu\min(u,v)}-1}{2\mu}\\
&=
\alpha^2\,\theta^2\,+\,\sigma^2\,\frac{e^{-\mu|u-v|}-e^{-\mu(u+v)}}{2\mu}\,.
\tag{10}
\end{align}
The last expression is known as the covariance function of an Ornstein-Uhlenbeck process.
Though some work is left to do it should be straightforward to perform the integral on the right hand side of (7) now.
