Showing that $X(t)=\int_0^tB(s)ds $ is a Gaussian process. The time integration of Brownian motion is given by $X(t)=\int_0^tB(s)ds$. In order to show that  $X(t)$ is a gaussian process I work with the definition that  for any $t_0<t_1 < t_2<\ldots <t_{n-1}$ the random variable
$$ \langle v,(X_{t_0}, ... , X_{t_{n-1}}) \rangle $$
is again Gaussian.
I tried to write out the quantity in terms of Riemann sums i.e.
$$ X(t) = \lim_{n\to {\infty}}\sum_{k=0}^{n-1}B(t_k)(t_{k+1}-t_k) $$
But got stuck with all the indices, how to continue?
We can assume that limits of gaussian random variables are again Gaussian.
 A: Using integration by parts we get:
$$X(t) = t B(t) - \int_0^ts dB(s)$$
Recall that the Itô integral of a deterministic function is a Gaussian process; thus $X(t)$ is the linear combination of two Gaussian process, making it Gaussian itself.

Edit: We can prove this without directly appealing to results from Itô calculus (as requested in the comments). Using summation by parts, we get: $$\sum_{k=0}^{n-1}B(t_{k+1})(t_{k+1} - t_k) = B(t_n)t_n - \sum_{k=1}^{n-1} t_k(B(t_{k+1})- B(t_{k}))$$
The right hand side can be seen to be a Gaussian process; sending $n \to \infty$ we get that $X(t)$ is a Gaussian process.
A: Maybe a more explicit way is sketched here, using directly the definition and w.o. integration by parts.
We have to show that $X(z_1),...,X(z_k)$ is jointly Gaussian, by definition of Gaussian process.
We can choose a partition of the time variable $u_{0}=0,...,u_{M}$ such that $u_{N_1}=z_1,...,u_{N_k}=z_k$ (that is, we have a partition that contains $z_1,..,z_k$). Call this partition $\Delta$.
The partition defines the approximations $X^{\Delta}(z_1),...,X^{\Delta}(z_k)$ in the standard way. We observe that these approximations are a finite combination of jointly distributed Gaussian ramdom variabless, therefore:

*

*$X^{\Delta}(z_1)=Z_0$ is Gaussian

*$X^{\Delta}(z_2)=X^{\Delta}(z_1)+Z_1$ , where $Z_1$ is Gaussian and $(Z_0,Z_1)$ have a joint Gaussian distribution

This shows that $(X^{\Delta}(z_1),X^{\Delta}(z_2))\sim(Z_0,Z_0+Z_1)$ is jointly Gaussian as well.
Analogously:

*

*$X^{\Delta}(z_3)=X^{\Delta}(z_2)+Z_2=X^{\Delta}(z_1)+Z_1+Z_2=$, where $Z_0,Z_1$,$Z_2$ are jointly Gaussian distributed.

This shows that $(X^{\Delta}(z_1),X^{\Delta}(z_2),X^{\Delta}(z_3)) \sim (Z_0,Z_0+Z_1,Z_0+Z_1+Z_2)$.
, and so on...
Now taking the limit $\Delta \rightarrow 0$ in $L^2$ (the one needed to define the Ito integral), if we know that the limit of a sequence of (multivariate) Gaussian variables is (multivariate) Gaussian (that we can take for granted according to the OP), we have shown that the resulting joint distribution of the limit process is Gaussian.
