Let $W(u)$ $(u \geq 0)$ be a Brownian motion on a probability space $(\Omega, \mathscr{F}, \mathbb{P})$. Let $I(T) = \int_0^T W(u) du$.
One can use integration by parts to show that $I(T) = \int_{0}^T (T-u) dW(u)$. Since the Ito integral of a deterministic integrand is normally distributed, one can use Ito's isometry to show that the variance of $I(T)$ is $\frac{1}{3}T^3$.
Page 112 of Monte Carlo Methods in Finance by Glasserman states that the variance of $I(T)$ is given by
$\mbox{Var}(I(T)) = 2\int_{0}^T \int_{0}^t \mbox{Cov}(W(u),W(t)) du dt = 2 \int_{0}^T \int_{0}^t u du dt = \frac{1}{3}T^3$.
I do not understand the first part of the equality - is it a consequence of Fubini's theorem?
Yes, it is a consequence of Fubini's theorem. By Fubini, we have
$$\text{var}(I(T)) = \mathbb{E} \left( \int_0^T W(u) \, du \int_0^T W_v \, dv \right) = \int_0^T \int_0^T \underbrace{\mathbb{E}(W_u W_v)}_{\text{cov}(W_u,W_v)} \, du \, dv.$$
Using again Fubini's theorem, we find
$$\begin{align*} \int_0^T \int_0^T \text{cov}(W_u,W_v) \, du \, dv &= \int_0^T \int_0^v\text{cov}(W_u,W_v) \, du \, dv + \int_0^T \int_v^T \text{cov}(W_u,W_v) \, du \, dv \\ &= \int_0^T \int_0^v\text{cov}(W_u,W_v) \, du \, dv + \int_0^T \int_0^u \text{cov}(W_u,W_v) \, dv \, du\\ &= 2 \int_0^T \int_0^r \text{cov}(W_r,W_t) \, dt \, dr \end{align*}$$
Combining both identities finishes the proof.