Let $e_1, e_2, e_3$ be a set of Cartesian basis vectors. Let $s = r e_1 + \theta e_2 + \phi e_3$ be a vector in some 3d vector space. If you can imagine, this is like having a stretched and deformed coordinate system where the spherical coordinates are instead rectangular.
This is clearly not the vector space we're used to, but there is a diffeomorphism--an invertible, differentiable map--between this vector space and the usual one. Let that map be denoted $f$:
$$f(s) = s' = r \sin \theta \cos \phi e_1 + r \sin \theta \sin \phi e_2 + r \cos \theta e_3 = x e_1 + y e_2 + z e_3$$
The Jacobian $\underline f$ of this map (which I will not calculate) gives the coordinate basis vectors of Euclidean space with respect to spherical coordinates.
$$\begin{align*} \underline f(e_1) &= e_r = \hat r \\ \underline f(e_2) &= e_\theta = r \hat \theta \\ \underline f(e_3) &= e_\phi =( r \sin \theta) \hat \phi \end{align*}$$
When you write something like $d\ell$ in an integral, you're writing it in terms of these spherical basis vectors $e_r, e_\theta, e_\phi$. You are not writing them in terms of unit vectors. Rather, you have
$$d\ell = e_r \, dr + e_\theta \, d\theta + e_\phi \, d\phi$$
If $\ell(t)$ is some curve, then what you're saying is that
$$\frac{d\ell}{dt} = e_r \frac{dr}{dt} + e_\theta \frac{d\theta}{dt} + e_\phi \frac{d\phi}{dt}$$
The quantity $d\ell/dt$ is a "tangent vector". It is literally the vector tangent to the curve at a given point. This is the key to understand about integration in multiple variables: we're always using some curve or surface or whatever and integrating by considering a small piece that is tangent to the curve or surface or whatever and summing over many such small pieces.
Volume integrals are no different. At each point, we take an infinitesimal volume. We build up these volumes by using three curves through a point, so their tangent vectors form a basis for the 3d space. The most convenient thing to do is to use the coordinate lines themselves.
Finally, we build up a representation of a tangent volume using what's called a "wedge product". That is, we define
$$dV \equiv (e_r \, dr) \wedge (e_\theta \, d\theta) \wedge (e_\phi \, d\phi) = (e_r \wedge e_\theta \wedge e_\phi) dr \, d\theta \, d\phi$$
(Note: in differential forms, the quantity $e_r \, dr$ is sometimes confusingly called $dr$ by itself. This, to me, is massively misleading.)
Wedge products are made to build up generalizations of vectors called $k$-vectors. The 3-vector $e_r \wedge e_\theta \wedge e_\phi$ can be written in terms of the unit vectors instead:
$$e_r \wedge e_\theta \wedge e_\phi = r^2 \sin \theta (\hat r \wedge \hat \theta \wedge \hat \phi) = \epsilon r^2 \sin \theta$$
$\epsilon$ is seen as the unit right-handed 3-vector. This choice of symbol is intentional; the components of $\epsilon$ are those of the Levi-Civita tensor. Any three unit vectors that are orthogonal form through wedges a 3-vector that is either $\pm \epsilon$. In particular, if you used Cartesian coordinates instead, you get
$$dV = (e_x \, dx ) \wedge (e_y \, dy) \wedge (e_z \, dz) = (\hat x \wedge \hat y \wedge \hat z) \, dx \, dy \, dz = \epsilon (dx \, dy \, dz)$$
Clearly, then, we have
$$dV = \epsilon (dx \, dy \, dz) = \epsilon r^2 \sin \theta \, (dr \, d\theta \, d\phi)$$
In multivariable calculus, you're not even prepared to deal with an object like $\epsilon$, so it is often silently thrown away, and all you're offered in explanation is that the volume $dV$ transforms according to determinants of Jacobians and such and such. All very true, but by throwing away $\epsilon$, the geometric picture of what is being done--how a tangent volume is constructed using the coordinate basis vectors--is utterly lost.
Forming volumes with the coordinate basis vectors is the very origin of how volume elements transform in new coordinate systems. The wedge product is key to doing this; vector calculus merely imitates this using the scalar triple product.
Best Answer
Not exactly an answer to your question, but a bit long for a comment:
Since $r$ denotes the radius of your sphere, it's best to use a different symbol, such as $\rho$, to denote the "radius function" (i.e., the distance to the origin).
Be that as it may, you can avoid the nuisance of integrating the angular variables in your question because $\nabla \cdot F = 3\rho^{2}$ is constant on each spherical shell $S_{\rho}$. Consequently, $$ \iiint_{V} (\nabla \cdot F)\, dV = \int_{0}^{r} \left(\iint_{S_{\rho}} 3\rho^{2}\, dA\right) d\rho = \int_{0}^{r} 3\rho^{2} \left(\iint_{S_{\rho}} dA\right) d\rho = 3(4\pi)\int_{0}^{r} \rho^{4}\, d\rho. $$