I came across this thread while searching for a similar topic.
In Nualart's book (Introduction to Malliavin Calculus), it is asked to show that $\int_0^t B_s ds$ is Gaussian and it is asked to compute its mean and variance. This exercise should rely only on basic Brownian motion properties, in particular, no Itô calculus should be used (Itô calculus is introduced in the next chapter of the book).
Here's a proposal:
Using, as a simplification, the variable change $s=tu$, one has that $\int_0^t B_s ds=tU_t$ where $U_t=\int_0^1 B_{tu}du$. Using a Riemann sum, one can write:
$$
U_t=\lim_{n\to\infty}\frac{1}{n}\sum_{k=0}^nB_{t\frac{k}{n}}=\lim_{n\to\infty}\frac{1}{n}S_n
$$
Using a summation by parts, one can write $S_n$ as:
\begin{align*}
S_n&=nB_t -\sum_{k=0}^{n-1} k \left(B_{t\frac{k+1}{n}}-B_{t\frac{k}{n}}\right) \\
&=n\sum_{k=0}^{n-1}\left(B_{t\frac{k+1}{n}}-B_{t\frac{k}{n}}\right)-\sum_{k=0}^{n-1} k \left(B_{t\frac{k+1}{n}}-B_{t\frac{k}{n}}\right) \\
&= \sum_{k=0}^{n-1} (n-k) \left(B_{t\frac{k+1}{n}}-B_{t\frac{k}{n}}\right) \\
&= \sum_{k=0}^{n-1} (n-k)X_{n,k}
\end{align*}
where $X_{n,k} := B_{t\frac{k+1}{n}}-B_{t\frac{k}{n}}$
Using B.M properties, we have that $\mathrm{Var}(X_{n,k})=\frac{t}{n}$, and $X_{n,k}$ are independent (as B.M increments).
We then have:
\begin{align*}
\mathrm{Var}(\frac{1}{n}S_n)&=\frac{1}{n^2} \sum_{k=0}^{n-1} (k-n)^2 \mathrm{Var}(X_{n,k})\\
&= \frac{t}{n^3} \sum_{k=0}^{n-1} (n-k)^2 \\
&= \frac{t}{n^3} \sum_{k=1}^{n} k^2 \\
&= t\frac{n(n+1)(2n+1)}{6n^3} \\
&= \frac{t}{3} + o(\frac{1}{n})
\end{align*}
Since we have
$\mathrm{Var}(\int_0^t B_s ds)=t^2\mathrm{Var}(U_t)=t^2\lim_{n\to\infty}\mathrm{Var}(\frac{1}{n}S_n)$,
we can conclude that
$$
\mathrm{Var}(\int_0^t B_s ds)=\frac{t^3}{3}
$$