Shreve II page 151, the Cox-Ingersoll-Ross model is defined as $$dr_t=(\alpha-\beta r_t)dt+\sigma\sqrt{r_t}dW_t$$ By applying Ito's Lemma, we obtain \begin{align} r_t&=r_0e^{-\beta t}+\frac{\alpha}{\beta}(1-e^{-\beta t})+\sigma e^{-\beta t}\int_0^te^{\beta u}\sqrt{r_u}dW_u\\ &=\frac{\alpha}{\beta}+\Big(r_0-\frac{\alpha}{\beta}\Big)e^{-\beta t}+\sigma e^{-\beta t}\int_0^te^{\beta u}\sqrt{r_u}dW_u \end{align} Now for the variance of $r_t$, Shreve suggests that we set $$X_t=e^{\beta t }r_t$$ and apply Ito's Lemma to obtain $dX_t$, after which $d(X_t^2)$ may be found. Then $d(X_t^2)$ is integrated to obtain $X_t^2$, from which $r_t^2$ is found. Finally, the variance is derived from $$Var(r_t)=E(r_t^2)-(E(r_t))^2$$ My question is, why not take the variance of $r_t$ immediately, that is \begin{align} Var(r_t)&=\sigma^2 e^{-2\beta t}\int_0^te^{2\beta u}E(r_u)du\\ &=\sigma^2 e^{-2\beta t}\int_0^te^{2\beta u}\Big(\frac{\alpha}{\beta}+\Big(r_0-\frac{\alpha}{\beta}\Big)e^{-\beta u}\Big)du \end{align} and from here the integration is simple. It yields the same result as Shreve's method.
One possibility is that I have assumed $E(r_u(dW_u)^2)$=$E(r_u)E((dW_u)^2)$, implying independence between $r_t$ and $(dW_t)^2$.
Any help is appreciated.