Operator $\nabla$

My teacher told me to solve some physics problem where I need to find $\rho(r)$ by using: $$\frac{\rho(r).dV}{\varepsilon_0}=\vec E(r+dr)S(r+dr)-\vec E(r)S(r)$$ where $dV$ is a small hollow spherical volume, i.e. $\frac d{dr}\left({\frac43\pi r^3}\right)=4\pi r^2dr$ and $\vec E(r)$ is a known vector (like $\vec E(r)=ar^4\hat r$ and a is a constant, after solving which we get $\rho=6ar^3\varepsilon_0$), and $S(r)$ is the spherical surface area $4\pi r^2$.

I found another convinient formula in a good book: $$\vec\nabla.\vec E=\frac{\rho}{\varepsilon_0}$$

I know that: $$\vec\nabla=\sum\frac{\partial}{\partial x}\hat x$$ and in spherical coordinates: $$\vec\nabla=\frac{\partial}{\partial r}\hat r+\frac{\partial}{r\partial \theta}\hat \theta+\frac{\partial}{r\sin\theta\partial \phi}\hat \phi$$ So I tried: $$\vec\nabla.\vec E=\left(\frac{\partial}{\partial r}\hat r+\frac{\partial}{r\partial \theta}\hat \theta+\frac{\partial}{r\sin\theta\partial \phi}\hat \phi\right).(ar^4\hat r)=\frac{\partial}{\partial r}ar^4=4ar^3\ne6ar^3$$

After checking back I found: $$\vec\nabla.\vec E=\frac{\partial(r^2E_r)}{r^2\partial r}+...=6ar^3$$

How did $r^2$ come there in?

The divergence should be calculated taking into account that the unit vector themselves are functions of the coordinates! 