Let $Y_t = \int_0^t X_s \mathrm{d} s$.
The process $Z_t = (X_t, Y_t)$ is also an Ito process, with
$\mathrm{d} Z_t = X_t (-a, 1) \mathrm{d} t + (\sigma, 0) \mathrm{d} W_t$ and the initial condition $Z_0 = (x_0, 0)$.
The process $Z_t$ is Gaussian, you determined $X_t$ to be Gaussian, and $Y_t$ is Gaussian as a linear functional of $X_t$. Thus to determine the joint distribution at time $t$, one needs to compute $\mathbb{E}(Z_t)$ and $\mathbb{E}(Z_t \otimes Z_t )$.
This would be done using Ito's lemma. Let $\tilde{X}_t = \mathrm{e}^{a t} X_t$ and $\tilde{Z}_t = ( \tilde{X}_t, Y_t )$. Then
$\mathrm{d} \tilde{Z}_t = (0, X_t) \mathrm{d} t + ( \sigma \mathrm{e}^{a t}, 0 ) \mathrm{d} W_t$. That is $\mathbb{E}(\tilde{Z}_t) = \mathbb{E}(\tilde{Z}_0) + \int_0^t (0, \mathrm{E}(X_s)) \mathrm{d} s $.
That is the $\tilde{X}_t$ process is martingale, hence
$$
\mathbb{E}( \tilde{X}_t ) = \mathbb{E}( \tilde{X}_0 ) = \mathbb{E}( X_0 ) = x_0
$$
Thus, $\mathbb{E}(X_t) = x_0 \exp(-a t)$, and by Ito's lemma $\mathbb{E}(Y_t) = \int_0^t x_0 \exp(-a s) \mathrm{d} s = \frac{x_0}{a} \left(1 - \mathrm{e}^{-a t} \right) $.
By Ito's lemma, again
$$
\mathrm{d}( \tilde{X}_t^2, \tilde{X}_t Y_t, Y_t^2 ) = ( \mathrm{e}^{2 a t} \sigma^2, \mathrm{e}^{a t} X_t^2, 2 X_t Y_t ) \mathrm{d} t + \sigma ( \mathrm{e}^{2 a t} X_t, \mathrm{e}^{a t} Y_t, 0 ) \mathrm{d} W_t
$$
Thus
$$
( \mathbb{E}(\tilde{X}_t^2), \mathbb{E}(\tilde{X}_t Y_t), \mathbb{E}(Y_t^2) ) =
( x_0^2, 0, 0 ) + \int_0^t ( \sigma^2 \mathrm{e}^{2 a s}, \mathrm{e}^{a s} \mathbb{E}(X_s^2), 2 \mathbb{E}( X_s Y_s) ) \mathrm{d} s
$$
This results in
$$
\begin{eqnarray}
\mathbb{E}\left(X_t^2\right) &=& x_0^2 \mathrm{e}^{-2 a t} + \frac{\sigma^2}{2} \frac{1 - \mathrm{e}^{-2 a t}}{a} \\
\mathbb{E}\left(X_t Y_t\right) &=& x_0^2 \mathrm{e}^{-a t} \left( \frac{ 1-\exp(-a t)}{a} \right) + \frac{\sigma^2}{2} \left( \frac{ 1 - \exp(-a t)}{a} \right)^2 \\
\mathbb{E}\left(Y_t^2\right) &=& x_0^2 \left( \frac{ 1 - \exp(-a t)}{a} \right)^2 + \frac{\sigma^2}{2 a^3} \left( 2 a t - 4 \left( 1 - \mathrm{e}^{-a t} \right) + \left( 1 - \mathrm{e}^{-2 a t} \right) \right)
\end{eqnarray}
$$
Converting the second moment into central moments, we get that the joint distribution of $Z_t$ is the bivariate normal distribution with means
$$
(m_1(t), m_2(t)) = x_0 \left( \mathrm{e}^{-a t}, \frac{1- \mathrm{e}^{-a t}}{a} \right)
$$
and covariance matrix:
$$
\Sigma(t) = \sigma^2 \left(\begin{array}{cc}
\frac{1-\mathrm{e}^{-2 a t}}{2 a} & \frac{1}{2} \left( \frac{1-\mathrm{e}^{-a t}}{a} \right)^2 \\ \frac{1}{2} \left( \frac{1-\mathrm{e}^{-a t}}{a} \right)^2 &
\frac{4 e^{-a t}+2 a t-3-e^{-2 a t}}{2 a^3} \end{array} \right)
$$
Added: Notice, that as $a$ becomes smaller the covariance matrix approaches the result from
this related question:
$$
\Sigma(t) \sim \sigma^2 \left(
\begin{array}{cc}
t(1 - a t + o(a)) & \frac{1}{2}t^2 \left(1-a t + o(a)\right) \\
\frac{1}{2}t^2 \left(1-a t + o(a)\right) & t^3 \left( \frac{1}{3} - \frac{a t}{4} + o(a) \right)
\end{array}
\right)
$$
The expectation of the Ito integral $\mathbb{E}( \int_0^t \mathrm{e}^{a s} \mathrm{d} W_s )$ is zero as George already said.
To compute $\mathbb{E}( W_t \int_0^t \mathrm{e}^{a s} \mathrm{d} W_s )$, write $W_t = \int_0^t \mathrm{d} W_s$. Then use Ito isometry:
$$
\mathbb{E}( W_t \int_0^t \mathrm{e}^{a s} \mathrm{d} W_s ) = \mathbb{E}\left( \int_0^t \mathrm{d} W_s \cdot \int_0^t \mathrm{e}^{a s} \mathrm{d} W_s \right) = \int_0^t (1 \cdot \mathrm{e}^{a s}) \mathrm{d} s = \frac{\mathrm{e}^{a t} - 1}{a} \phantom{hhhh}
$$
Best Answer
Using polarization, Ito's isometry implies $ \mathbb{E}\left( \left( \int_0^t f(s) \mathrm{d} W_s \right) \left( \int_0^t g(s) \mathrm{d} W_s \right) \right) = \int_0^t f(s) g(s) \mathrm{d} s$.
Thus, for $s<t$ $$ \begin{eqnarray} \mathbb{E} \left( X_s X_t \right) &=& \mathbb{E} \left( \left( \sigma \int_0^t \mathrm{e}^{-a(s - \tau)} \mathbf{1}\left(\tau \le s\right) \mathrm{d}W_\tau \right) \left( \sigma \int_0^t \mathrm{e}^{-a(t - \tau)} \mathrm{d}W_\tau \right) \right) \\ &=& \sigma^2 \int_0^{\min(s,t)} \mathrm{e}^{-a (s+t-2 \tau)}\mathrm{d} \tau = \sigma^2 e^{-a (s+t)} \cdot \frac{e^{2 a \min (s,t)}-1}{2 a} = \sigma^2 e^{-a t} \cdot \frac{\sinh(a s)}{a} \end{eqnarray} $$
Using $W_s = \int_0^s 1 \,\, \mathrm{d} W_\tau$, we similarly get $$ \mathbb{E} \left( X_t W_s\right) = \sigma \int_0^s \mathrm{e}^{-a(t-\tau)}\mathrm{d} \tau = \sigma \mathrm{e}^{-a t} \cdot \frac{ \left(\mathrm{e}^{a s}-1\right)}{a} $$