First, since you're integrating the divergence of $E$, using $E$ itself would not make sense.
Second, as long as the divergence is not constant, well, you can't pull it out of the integral. You're basically suggesting you should be able to pull something out of the integral so that you can rewrite it as $(\nabla \cdot E) \int 1 \, dV$. You can only do this if the divergence is constant--if it doesn't depend on any of the integration variables $r, \phi, z$.
In short, how could you do this...
$\int_V \rho(r, \phi z) \, r \, dr \, d\phi \, dz \to \rho(r, \phi, z) \int_V r \, dr \, \phi \, dz$
...if integrating completely removes all dependence on $r, \phi, z$? You can't pull something out that depends on integration variables. That's the problem.
In your particular problem, however, you can use the divergence theorem to convert this volume integral into a surface integral. Then, you integrate $E$ directly on the surface of the cylinder. This can be convenient because then, for example, $E \cdot \hat n$ may be constant on the curved wall of the cylinder (this is often the case in problems involving cylindrical symmetry).