$\DeclareMathOperator{\div}{div}
\DeclareMathOperator{\tr}{tr}$
We consider the following definition of the divergence. Let $(M,g)$ be an orientable Riemannian manifold, with Riemannian volume form $\mathrm{d}vol_g$, and $X$ be a vector field. The divergence of $X$ is the unique smooth function $\div X$ such that $L_X\mathrm{d}vol_g = (\div X) \mathrm{d}vol_g$.
Let $\{e_1,\ldots,e_n\}$ be a local orthonormal frame and consider its dual frame $\{\theta^1,\ldots,\theta^n\}$, that is $\theta^i = g(\cdot,e_i)$. Then we have the equality
$$
\mathrm{d}vol_g = \theta^1\wedge\cdots\wedge \theta^n.
$$
Indeed, these two volume forms are proportional and coincide on the orthonormal frame $\{e_1,\ldots,e_n\}$. Now, since $L_X$ is a derivation of the exterior algebra, it follows that
$$
L_X\mathrm{d}vol_g = \sum_{i=1}^n \theta^1\wedge\cdots\wedge L_X\theta^i \wedge \cdots \wedge \theta^n.
$$
For any vector field $Y$, we have
\begin{align}
(L_X\theta^i)(Y) &= X\theta^i(Y) - \theta^i([X,Y]) & \text{by Leibniz rule}\\
&= Xg(Y,e_i) - g([X,Y],e_i)\\
&= g(\nabla_XY-[X,Y],e_i) + g(Y,\nabla_Xe_i) &\text{from the compatibility of $g$ and $\nabla$}\\
&= g(\nabla_YX,e_i) + g(Y,\nabla_Xe_i) & \text{since $\nabla$ is torsion-free}\\
&= (\theta^i\circ \nabla X)(Y) + g(Y,\nabla_Xe_i).
\end{align}
Applying $\nabla_X$ to the equality $\|e_i\|^2 = 1$ gives $\nabla_Xe_i\perp e_i$, so that the 1-form $g(\cdot,\nabla_Xe_i)$ vanishes on $e_i$, and is then a linear combination of $\{\theta^j\}_{j\neq i}$. It then disappears in the wedge product, and it follows that
$$
L_X\mathrm{d}vol_g = \sum_{i=1}^n \theta^1\wedge\cdots\wedge (\theta^i\circ \nabla X) \wedge \cdots \wedge \theta^n.
$$
Finally, evaluating on the orthonormal frame $\{e_1,\ldots,e_n\}$ yields
\begin{align}
\div X &= \sum_{i=1}^n \theta^1(e_1)\times \cdots \times \theta^i(\nabla_{e_i}X) \times \cdots \theta^n(e_n)\\
&= \sum_{i=1}^n \theta^i(\nabla_{e_i}X)\\
&= \sum_{i=1}^n g(\nabla_{e_i}X,e_i)\\
&= \tr \nabla X.
\end{align}
There are plenty of other proofs (on this website also) relying on the coordinate expression of the Riemannian volume form. This method is known as the method of moving frames (referring to the local orthonormal (co-)frame). It surprisingly turns out that it is less known whilst it is (in my opinion) simpler. Some people (including me) prefer to work in a coordinate-free manner, and to that end, moving frames are good friends.
As it has been mentioned in comment, by Bochner's formula for function $f$ the answer is clear:
$$\Delta\nabla_if =
\nabla_i\Delta f + \mbox{R}_{ij}\nabla_jf.$$
Obviously in flat metrics we have $\mbox{R}_{ij}=0$ and ..
Best Answer
This calculation is done in full detail (without the assumption that the basis is orthonormal) in example $10.1.32$ of Liviu Nicolaescu's wonderful book Lectures on the Geometry of Manifolds which is freely available as a pdf from his website.