$\DeclareMathOperator\div{div}$The formula for $\nabla\cdot X$ is incorrect. The notation with the 'usual' dot product is misleading. Properly it is:
$$\div F = \frac 1\rho\frac{\partial(\rho F^i)}{\partial x^i}$$
where $\rho=\sqrt{\det g}$ is the coefficient of the differential volume element $dV=\rho\, dx^1\wedge\ldots \wedge dx^n$, meaning $\rho$ is also the Jacobian determinant, and where $F^i$ are the components of $F$ with respect to an unnormalized basis.
In polar coordinates we have $\rho=\sqrt{\det g}=r$, and:
$$\div X = \frac 1r \frac{\partial(r X^r)}{\partial r}
+ \frac 1r\frac{\partial(r X^\theta)}{\partial \theta}$$
In the usual normalized coordinates $X=\hat X^{r}\frac{\partial}{\partial r} + \hat X^{\theta}\frac 1r\frac{\partial}{\partial\theta}$ this becomes:
$$\div X = \frac 1r \frac{\partial(r \hat X^{r})}{\partial r}
+ \frac 1r\frac{\partial \hat X^{\theta}}{\partial \theta}$$
which agrees with the usual formula given in calculus books.
First, you have to understand that not all bases are coordinate bases. A coordinate basis is one where the basis vectors are of the form $e_\mu = \partial_\mu = \frac{\partial}{\partial x^\mu}$ for some coordinate system $\{x^\mu\}_{\mu = 1}^d$. Equivalently, a coordinate basis is one in which the commutator $[e_\mu, e_\nu]$ between any pair of basis vectors vanishes. In differential geometry, ON bases (in which the metric is $g_{\mu\nu} = \delta_{\mu\nu}$) have a tendency not to be coordinate bases!
The regular polar coordinate system $(r,\theta)$ gives you a basis $\{\partial_r, \partial_\theta\}$ and a metric $g = \begin{pmatrix} 1 & 0 \\ 0 & r^2\end{pmatrix}$ with a determinant $\sqrt{\det(g)} = r$. In this basis, you would express a vector as $V = V^\mu \partial_\mu = V^r \partial_r + V^\theta \partial_\theta$.
Since this is a coordinate basis, we can use the given formula to evaluate the divergence:
$$ \nabla_aV^a = \nabla_\mu V^\mu = \frac{1}{\sqrt{\det(g)}} \partial_\mu \left(\sqrt{\det(g)} V^\mu\right) = \frac{1}{r} \partial_\mu (rV^\mu) = \underline{\partial_r V^r + \frac{1}{r} V^r + \partial_\theta V^\theta}. $$
If you switch your basis to $\{e_1 = \partial_r, e_2 = \frac{1}{r} \partial_\theta\}$, the metric becomes $g_{\mu\nu} = \delta_{\mu\nu}$ and a vector is expressed as $V = V^\mu e_\mu = V^1 e_1 + V^2 e_2$, from which you can see that the new vector components are $V^1 = V^r$ and $V^2 = r V^\theta$. The new basis is not a coordinate basis; you can easily check that $[e_1, e_2] = -\frac{1}{r^2} \partial_\theta \neq 0$.
The expression $(\nabla_\mu V)^\nu = \partial_\mu V^\nu + \Gamma^\nu_{\rho\mu} V^\rho$ only holds in a coordinate basis. The symbol $\partial_\mu$ has no meaning in a non-coordinate basis. Instead, the general formula is
$$ (\nabla_\mu V)^\nu = e_\mu(V^\nu) + \Gamma^\nu_{\rho\mu} V^\rho. $$
Here, $e_\mu(V^\nu)$ means that the differential operator $e_\mu$ acts on the function $V^\nu$. We cannot use the given formula for the divergence in this basis, because it only makes sense in a coordinate basis. Furthermore, even though the metric is $\delta_{\mu\nu}$ in this basis, the connection components $\Gamma^\mu_{\nu\rho}$ do not vanish since their expression as Christoffel symbols $\Gamma^\mu_{\nu\rho} = \frac{1}{2} g^{\mu\sigma} (g_{\nu\sigma,\rho} + g_{\sigma\rho,\nu} - g_{\nu\rho,\sigma})$ only makes sense in a coordinate basis.
To find the connection components in a general basis, we must use the general expression for the Levi-Civita connection,
$$ Z \cdot \nabla_X Y = \frac{1}{2}\Big(X(Y \cdot Z) + Y(X \cdot Z) - Z(X \cdot Y) - X \cdot [Y,Z] - Y \cdot [X,Z] + Z \cdot [X,Y]\Big). $$
The connection components are defined as $\Gamma^\mu_{\nu\rho} = (\nabla_\rho e_\nu)^\mu = g^{\mu\sigma} (\nabla_\rho e_\nu)_\mu = g^{\mu\sigma} e_\mu \cdot \nabla_{e_\rho} e_\nu$, so we find (with $X = e_\rho$, $Y = e_\nu$, $Z = e_\sigma$)
$$ \Gamma^\mu_{\nu\rho} = \frac{1}{2} g^{\mu\sigma} \Big(
e_\rho (\overbrace{e_\nu \cdot e_\sigma}^{g_{\nu\sigma}})
+ e_\nu (\overbrace{e_\rho \cdot e_\sigma}^{g_{\rho\sigma}})
- e_\sigma (\overbrace{e_\nu \cdot e_\rho}^{g_{\nu\rho}})
- e_\rho \cdot [e_\nu, e_\sigma]
- e_\nu \cdot [e_\rho, e_\sigma]
+ e_\sigma \cdot [e_\nu, e_\rho]
\Big). $$
Let us use this formula to compute $\nabla_\mu V^\mu$ in the basis $\{e_1, e_2\}$. Since the metric components in our case are constant, the first three terms vanish (these are the terms that give the Christoffel symbols in a coordinate basis). But the other three terms do not vanish. Inserting $g^{\mu\sigma} = \delta^{\mu\sigma}$, we get
$$ \Gamma^\mu_{\nu\rho} = -\frac{1}{2} \Big(
e_\rho \cdot [e_\nu, e_\mu]
+ e_\nu \cdot [e_\rho, e_\mu]
- e_\mu \cdot [e_\nu, e_\rho]
\Big). $$
As we found before, we have the commutator $[e_1, e_2] = -\frac{1}{r^2} \partial_\theta = -\frac{1}{r} e_2$. To compute $\nabla_\mu V^\mu = e_\mu(V^\mu) + \Gamma^\mu_{\nu\mu} V^\nu$, we need $\Gamma^\mu_{1\mu}$ and $\Gamma^\mu_{2\mu}$. Clearly $\Gamma^\mu_{\mu\mu}$ (no sum) vanishes because all commutators vanish, so $\Gamma^\mu_{1\mu} = \Gamma^2_{12}$ and $\Gamma^\mu_{2\mu} = \Gamma^1_{21}$. Doing the computation:
$$ \begin{align}
\Gamma^2_{12} &= -\frac{1}{2} \Big(e_2 \cdot [e_1, e_2] + e_1 \cdot [e_2, e_2] - e_2 \cdot [e_2, e_1]\Big) = -e_2 \cdot [e_1, e_2] = \frac{1}{r} e_2 \cdot e_2 = \frac{1}{r}, \\
\Gamma^1_{21} &= -\frac{1}{2} \Big(e_1 \cdot [e_2, e_1] + e_2 \cdot [e_1, e_1] - e_1 \cdot [e_1, e_2]\Big) = -e_1 \cdot [e_2, e_1] = -\frac{1}{r} e_1 \cdot e_2 = 0.
\end{align} $$
We are then able to calculate the divergence:
$$ \nabla_\mu V^\mu = e_1(V^1) + \overbrace{\Gamma^1_{21}}^0 V^2 + e_2(V^2) + \overbrace{\Gamma^2_{12}}^{1/r} V^1 $$
and inserting the expressions for $e_1$, $e_2$, $V^1$ and $V^2$:
$$ \begin{align}
\nabla_\mu V^\mu &= \partial_r(V^1) + \frac{1}{r} V^1 + \frac{1}{r} \partial_\theta(V^2) \\
&= \partial_r (V^r) + \frac{1}{r} V^r + \frac{1}{r} \partial_\theta (r V^\theta) \\
&= \underline{\partial_r V^r + \frac{1}{r} V^r + \partial_\theta V^\theta},
\end{align} $$
which is the same as before.
It is, in fact, always possible to find an orthonormal coordinate basis, at least locally. This is what is called a local inertial frame in GR: Starting from an ON basis $\{e_\mu\}_{\mu = 1}^d$ at a point $p$, you can define coordinates (normal coordinates) around $p$ via the geodesics through $p$. Starting from the basis $\{\partial_r, \partial_\theta\}$, we would get the regular Cartesian coordinate system, but centred on our starting point $p$, given by $(r_0, \theta_0)$:
$$ \begin{align}
\tilde{x} &= (r - r_0) \cos(\theta - \theta_0) \\
\tilde{y} &= (r - r_0) \sin(\theta - \theta_0).
\end{align} $$
In this coordinate system, you will indeed find a diagonal metric as well as vanishing Christoffel symbols (generally only true locally, but here globally), which means that
$\nabla_\mu V^\mu = \partial_\mu V^ \mu$.
Best Answer
The polar coordinates orthonormal basis is $\partial_\mu = \{\partial_r, \frac{1}{r}\partial_\theta\}$, which you can see by transforming from the cartesian coordinate system. Hence the $\frac{1}{r}$ in the final expression.