Prologue:
My definition of standard Brownian motion is a (nonstationary) Gaussian random process $\{W_t \colon t > 0\}$ where
$$W_t = \int_0^t X_s \,\mathrm ds,$$
$\{X_s \colon s \in \mathbb R \}$ being a (zero-mean) white Gaussian noise process with autocorrelation function $\sigma^2 \delta(s)$ where $\delta(s)$ denotes the Dirac delta.Thus, we have that
$$E[W_t] = E\left[\int_0^t X_s \,\mathrm ds \right] = \int_0^t E[X_s] \,\mathrm ds = 0,$$
and
\begin{align}
E[(W_t)^2] &= E\left[\int_0^t X_s \,\mathrm ds \int_0^t X_r \,\mathrm dr \right]\\
&= E\left[\int_0^t \int_0^t X_sX_r \,\mathrm ds\,\mathrm dr \right]\\
&= \int_0^t \int_0^t E[X_sX_r] \,\mathrm ds\,\mathrm dr\\
&= \int_0^t \int_0^t \sigma^2\delta(s-r) \,\mathrm ds\,\mathrm dr\\
&= \int_0^t \sigma^2\,\mathrm dr\\
&= \sigma^2 t.
\end{align}
To find $E[W_tW_s]$, note that
$$W_tW_s = W_{\min(t,s)}W_{\max(t,s)} = W_{\min(t,s)}\left(W_{\min(t,s)}+\int_{\min(t,s)}^{\max(t,s)}X_r \, \mathrm dr\right)$$
where $W_{\min(t,s)}$ and $\displaystyle\int_{\min(t,s)}^{\max(t,s)}X_r \, \mathrm dr$, being the integrals of white noise over two disjoint intervals of time, thus are independent zero-mean random variables. Consequently,
\begin{align} E[W_tW_s] &= E\left[W_{\min(t,s)}\left(W_{\min(t,s)}+\int_{\min(t,s)}^{\max(t,s)}X_r \, \mathrm dr\right)\right]\\
&= E[(W_{\min(t,s)})^2] + E\left[W_{\min(t,s)}\int_{\min(t,s)}^{\max(t,s)}X_r \, \mathrm dr\right)]\\
&= E[(W_{\min(t,s)})^2] = \sigma^2\min(t,s)
\end{align}
Summary: Standard Brownian motion is a zero-mean nonstationary Gaussian process with autocorrelation function $R_W(t,s) = E[W_tW_s] = \sigma^2\min(t,s).$
End of Prologue:
Since $E[W_tW_s]=\sigma^2\min(t,s)$,
\begin{align}
E\left[\left(\int_0^T W_t\,\mathrm dt\right)^2\right]&=E\left[\int_0^TW_t\,\mathrm dt\int_0^TW_t\,\mathrm dt\right]\\
&=E\left[\int_0^TW_t\,\mathrm dt\int_0^TW_s\,\mathrm ds\right]\\
&=E\left[\int_0^T\int_0^TW_t\,W_s\,\mathrm dt \,\mathrm ds\right]\\
&=\int_0^T\int_0^T E[W_tW_s]\,\mathrm dt\,\mathrm ds\\
&=\int_0^T\int_0^T \sigma^2\min(t,s)\,\mathrm dt\,\mathrm ds\\
&=\int_0^T\left[\int_0^T\sigma^2\min(t,s)\,\mathrm dt\right]\,\mathrm ds\\
&=\int_0^T\left[\int_0^s\sigma^2 \min(t,s)\,\mathrm dt
+ \int_s^T\sigma^2 \min(t,s)\,\mathrm dt\right]\,\mathrm ds\\
&=\int_0^T \int_0^s\sigma^2 \min(t,s)\,\mathrm dt\,\mathrm ds
+ \int_0^T\int_s^T\sigma^2 \min(t,s)\,\mathrm dt\,\mathrm ds\tag{1}\\
&= \int_0^T\int_0^s\sigma^2 t\,\mathrm dt\,\mathrm ds
+ \int_0^T\int_s^T\sigma^2 s\,\mathrm dt\,\mathrm ds\tag{2}\\
&=\int_0^T\sigma^2\frac{s^2}{2}\,\mathrm ds + \int_0^T\sigma^2(T-s)s\,\mathrm ds\\
&= \sigma^2\left(\frac{T^3}{6}+\frac{T^3}{2}-\frac{T^3}{3}\right)\\
&= \sigma^2\frac{T^3}{3}<\infty
\end{align}
I hope that the answer above is understandable and acceptable to most readers, but for the benefit of @youpilat13 and any others confused as to how $(2)$ follows from $(1)$, I add the following explanation.
There are two double integrals in $(1)$ and they are over disjoint right triangular regions, separated by the line $s=t$, in the plane. As to why the two integrals are being summed, we are trying to compute the integral over a square region and we are breaking up the integral over the square into the sum of integrals over two disjoint subsets (right triangles), computing each integral separately, and then adding the results together to determine the integral over the square region.
In the inner integral in the first double integral in $(1)$, $s$ is a constant with respect to $t$, the variable of integration, but since $t$ varies from $0$ to $s$ only, it must be that $\min(t,s)$ equals $t$. Similarly, in the inner integral in the second double integral in $(1)$, $s$ is a constant with respect to $t$, the variable of integration, but since $t$ varies from $s$ to $T$ only, it must be that $\min(t,s)$ equals $s$.
