How to show that the divergence can be written as the trace of total covariant derivative, ${\rm div}={\rm tr}\circ\nabla$?

We can show that this equality holds at any point of $M$ by computing in local coordinates. At the point $p\in M$ we use Riemannian normal coordinates, which have the property that the metric at $p$ is given by the identity matrix, $g_{ij}(p)=\delta_{ij}$, and that the derivatives of the metric coefficients are zero at $p$, as are the Christoffel symbols. With this we can write \begin{equation} \text{tr}\ \nabla X (p) = \sum_{i=1}^n g(\nabla_{\partial_i} X,\partial_i)(p) = \sum_{i=1}^n(\partial_i X^j)(p)\delta_{ij} = \sum_{i=1}^n\partial_i X^i(p). \end{equation} On the other hand, we have $\text{det}(g)(p)=1$ in these coordinates, and since the determinant is polynomial in the entries of the matrix of $g$, all its partial derivatives vanish at $p$ in these coordinates. With $dV_g = \sqrt{|g|}dx^1\wedge\ldots\wedge dx^n$ we can write \begin{equation} X\lrcorner\ dV_g = \sqrt{|g|}\sum_{i=1}^n(-1)^{i+1}X^i dx^1\wedge\ldots\widehat{dx^i}\wedge\ldots\wedge dx^n, \end{equation} where the hatted factor is the one excluded from the wedge product. Noting that the partial derivatives of $\sqrt{|g|}$ vanish at $p$ we can write \begin{equation} d(X\lrcorner\ dV_g)(p) = \sum_{i=1}^n(-1)^{i+1}(\partial_j X^i)(p)dx^j\wedge dx^1\wedge\ldots\widehat{dx^i}\wedge\ldots\wedge dx^n. \end{equation} All the terms where $i\neq j$ will be zero by antisymmetry of $\wedge$. Permuting the leftmost factor $dx^j$ to the place of the hatted factor we change sign $i-1$ times, so the sign factors cancel. We are left with \begin{equation} d(X\lrcorner\ dV_g)(p) = \left(\sum_{i=1}^n(\partial_i X^i)(p)\right)\ dx^1\wedge\ldots\wedge dx^n = \text{tr}\ \nabla X(p)\ dV_g. \end{equation}

I'll add another completely different solution. If $M$ is any manifold, $\nabla$ is any torsion-free connection, $X$ is any vector field and $\beta\in \Omega^k(M)$, then $$\mathscr{L}_X\beta = \nabla_X\beta +\sum_{i=1}^k\beta(\cdot,\ldots,\nabla X,\ldots,\cdot).$$So if $\omega\in\Omega^n(M)$ is a top degree form, we have $$\mathscr{L}_X\omega = \nabla_X\omega+({\rm div}\,X)\omega$$by linear algebra. Cartan's magic formula turns that into $${\rm d}(\iota_X\omega) = \nabla_X\omega+({\rm div}\,X)\omega.$$If $\nabla$ is the Levi-Civita connection of a pseudo-Riemannian manifold $(M,g)$ and we apply this to ${\rm d}V_g$, which is parallel (as the Ricci tensor of the Levi-Civita connection is symmetric), it follows that $${\rm d}(\iota_X{\rm d}V_g)=({\rm div}\,X){\rm d}V_g$$as wanted. So this is a more powerful proof that does not rely on local computations, and holds for metrics of any signature.

We can work locally since everything there is a tensor. Assume $(E_1,\ldots,E_n)$ is a local orthonormal frame field, such that $(\nabla_{E_j}E_i)_p = 0$ for some $p \in M$. We will check the equality at $p$.

If $(\theta^1,\ldots, \theta^n)$ is the corresponding coframe, then $$dV_g = \theta^1 \wedge \cdots \wedge \theta^n.$$We can write $X = \sum_{i=1}^n \theta^i(X) E_i$, and we also know that $\theta^i(X) = \langle X,E_i\rangle$, by orthonormality. On one hand, it is easy (?) to check that $${\rm tr}(\nabla X) = \sum_{i=1}^n \langle \nabla_{E_i}X,E_i\rangle,$$so let's compute the other expression. We start computing the interior derivative $$\begin{align}\iota_X(dV_g) &= \iota_X(\theta^1 \wedge \cdots \wedge \theta^n) \\ &= \sum_{i=1}^n (-1)^{i-1} \theta^1 \wedge \cdots \wedge \iota_X(\theta^i)\wedge \cdots \wedge \theta^n \\ &= \sum_{i=1}^n (-1)^{i-1} \theta^i(X) \;\theta^1\wedge \cdots \wedge\widehat{\theta^i} \wedge \cdots \wedge \theta^n, \end{align}$$where the hat denotes omission. Then we compute the exterior derivative of that: $$\begin{align} d(\iota_X(dV_g)) &= \sum_{j=1}^n \sum_{i=1}^n (-1)^{i-1} E_j(\theta^i(X)) \; \theta^j \wedge \theta^1\wedge \cdots \wedge\widehat{\theta^i} \wedge \cdots \wedge \theta^n \\ &= \sum_{i=1}^n (-1)^{i-1} E_i(\theta^i(X)) \;\theta^i\wedge \theta^1\wedge \cdots \wedge\widehat{\theta^i} \wedge \cdots \wedge \theta^n \\ &= \sum_{i=1}^n E_i(\theta^i(X)) \;\theta^1\wedge \cdots \wedge \theta^n \\ &= \left(\sum_{i=1}^n E_i(\theta^i(X))\right) \theta^1\wedge \cdots \wedge \theta^n \\ &= \left(\sum_{i=1}^n E_i(\theta^i(X))\right) dV_g.\end{align}$$The result follows since at $p$ we have $$(E_i)_p(\theta^i(X)) = \langle (\nabla_{E_i}X)_p,(E_i)_p\rangle + \langle X, (\nabla_{E_i}E_i)_p\rangle =\langle (\nabla_{E_i}X)_p,(E_i)_p\rangle. $$The result is also true for pseudo-Riemannian manifolds, with a slight adaptation of the calculations done above.