I came across this thread while searching for a similar topic.
In Nualart's book (Intoduction to Malliavin Calculus), it is asked to show that $\int_0^t B_s ds$ is gaussien and it is asked to compute its mean and variance. This exerice should rely only on basic brownian motion properties, in particular, no Itô calculus should be used (Itô calculus is introduced in the next cahpter 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*}
nS_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 propertries, we have that $\mathrm{Var}(X_{n,k})=\frac{t}{n}$, and $X_{n,k}$ are independant (as B.M increments).
We then have:
\begin{align*}
\mathrm{Var}(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)$, we can conclude that
$$
\mathrm{Var}(\int_0^t B_s ds)=\frac{t^3}{3}
$$