I am eager to compute the following variance, where $W_t$ is a Wiener process:
$$\mathbb{V}\left[I_t\right], \text{ where } I_t = \int\limits_0^tW_u^2 u^2\partial u.$$
I have already computed $\mathbb{E}[I_t] = \Large\frac{t^4}{4}$, and hence currently I have:
$$ \mathbb{V}[I_t] = \mathbb{E}[I^2_t] - \mathbb{E}^2[I_t] = \mathbb{E}[I^2_t] - \frac{t^{16}}{16} = \mathbb{E}\left[\left(\int\limits_0^tW_u^2 u^2\partial u\right)^2\right] - \frac{t^{8}}{16}. $$
And here, I do not understand, how to compute the first summand. I thought of applying the Itô isometry here, but as long as $W_u^2u^2$ is a stochastic process (at least it seems for me it is, what may be verified), the differential is taken not over a Wiener process, but rather over the time-like variable $u$, so that I have no further ideas on how to proceed.
Any help would be appreciated, thank you in advance!