Using the definition of a Riemann integral, we have that the Riemann sum for a partition $0=t_0<t_1<\dots<t_N=T$ is
$$
S_N = \sum_{i=1}^N W(t_i)(t_i - t_{i-1}),
$$
Now, we have that $W(t_i) = \sum_{k=1}^i[W(t_k) - W(t_{k-1})]$, so
$$
S_N = \sum_{i=1}^N\sum_{k=1}^i[W(t_k) - W(t_{k-1})](t_i - t_{i-1}).
$$
$$
S_N = \sum_{i=1}^N[W(t_i) - W(t_{i-1})](T - t_{i-1}).
$$
Since the increments of the Wiener process are independent and Gaussian with variance equal to the time increment, we have that
$$
E[S_N^2] = \sum_{i=1}^N[t_i - t_{i-1}](T - t_{i-1})^2.
$$
The expression above is a Riemann sum too, so when the mesh of the partition vanishes and $N\to\infty$, we have that
$$
E\left[\left(\int_0^T W(t)\,dt\right)^2\right] = \int_0^T (T-t)^2\,dt = \frac{T^3}{3}
$$
Gaussian random variables form a linear subspace of $L_2$, since Gaussianity is closed with respect to addition and multiplication by scalar, i.e., if X and Y are Gaussian random variables then
$$
Z = aX + bY
$$
is also a Gaussian random variable. A subspace of a Hilbert space is also a Hilbert space, so it is complete. You can then use dominated convergence to show that $S_N$ converges in $L_2$ (mean square convergence), so that its limit is also Gaussian. Taking the expected value of $S_N$ it is also easy to see that the integral is zero-mean.