Let's be a bit systematic. Let $M$ be a smooth $m$-dimensional manifold. Suppose that $M$ is orientable and oriented (so that we can use $m$-forms as densities). Let $\xi:E\rightarrow M$ be a smooth vector bundle. The adjoint bundle of $(E,\xi,M)$ is denoted $(E^+,\xi^+,M)$ with $$ E^+=E^\ast\otimes\Lambda^m(M), $$ and we identify $E^+_x\cong\mathrm{Hom}(E_x,\Lambda^m_x(M))$ and $\Gamma(E^+)\cong\mathrm{Hom}_{C^\infty(M)}(\Gamma(E),\Omega^m(M))$, i.e. an element of a fibre of $E^+$ is interpreted as a linear map from the fibre $E_x$ to the fibre $\Lambda^m_x(M)$ of $m$-forms at $x$, and a smooth section of $E^+$ is interpreted as a $C^\infty(M)$-linear map from sections of $E$ to $m$-forms.
Given a pair $\xi,\zeta:E,F\rightarrow M$ of vector bundles, and a linear differential operator $D:E\rightarrow F$ (note that this is an abuse of notation because a linear diff. op. is not actually a map between the total spaces of vector bundles), the formal adjoint of $D$ is a linear differential operator $D^+:F^+\rightarrow E^+$ determined uniquely by the relation $$ \langle\psi,D\phi\rangle-\langle D^+\psi,\phi\rangle = \mathrm{d}\Xi(\phi,\psi), $$ where $\phi\in\Gamma(E)$, $\psi\in\Gamma(F)$, $\langle -,-\rangle$ is the $m$-form valued (!!) pairing of a vector bundle and its adjoint bundle, and $\Xi:E\times F\rightarrow \Lambda^{m-1}(M)$ is a bivariate linear differential operator.
Suppose that $U\subseteq M$ is a coordinate chart domain in $M$ which is also a local frame domain for both $E$ and $F$, and the linear differential operator takes the coordinate form $$ (D\phi)^\alpha=\sum_{k=0}^rD^{\alpha,i_1...i_k}_a\partial_{i_1}\dots\partial_{i_k}\phi^a, $$ then the formal adjoint has coordinate form $$ (D^+\psi)_a=\sum_{k=0}^r(-1)^k\partial_{i_1}\dots\partial_{i_k}\left(D^{\alpha,i_1...i_k}_a\psi_\alpha\right). $$
Given a vector bundle $\xi:E\rightarrow M$ the $\mathbb R$-linear space $\mathcal D^\sharp(E)$ of distributions with values in $E$ is defined as the set of all (appropriately continuous) linear functions $\mathcal D(E)\rightarrow\mathbb R$, where $\mathcal D(E):=\Gamma_c(E^+)$ is the set of all smooth, compactly supported sections of the adjoint bundle $E^+$.
Then if $D:E\rightarrow F$ is any linear differential operator and $T\in\mathcal D^\sharp(E)$ is a distribution, the linear differential operator acts on $T$ through the formal adjoint, i.e. $$ (DT)(\varphi):=T(D^+\varphi),\quad\forall\varphi\in \mathcal D(F). $$
The Dirac delta $\delta_p\in\mathcal D^\sharp(\Lambda^m(M))$ is a distributional $m$-form, defined on any test function $\varphi\in\mathcal D(\Lambda^m(M))\cong C^\infty_c(M)$ as $$ \delta_p(\varphi):=\varphi(p). $$
The covariant derivative $\nabla$ associated to any smooth linear connection is a linear differential operator of order $\le 1$ on any tensor bundle. In our particular case, it is a differential operator $\nabla:\Lambda^m(M)\rightarrow T^\ast(M)\otimes\Lambda^m(M)\cong T^+(M),$ so its formal adjoint is of the type $\nabla^+:T(M)\rightarrow F(M)$, where $F(M)=M\times\mathbb R$ is the trivial line bundle.
For any $m$-form $\omega\in \Omega^m(M)$ with $\omega=\rho\,\mathrm d^mx$, we have $$ \nabla\omega=\left(\partial_i\rho-\Gamma^{j}_{ji}\rho\right)\mathrm{d}x^i\otimes\mathrm{d}^mx, $$ so applying the formula for the formal adjoint, we get $$ \nabla^+X=-\partial_i X^i-\Gamma^{j}_{ji}X^i. $$
So then, we have $$ \nabla\delta_p (X)=\delta_p(\nabla^+X)=(\nabla^+X)(p)=-\partial_i X^i(p)-\Gamma^{j}_{ji}(p)X^i(p).$$
On the other hand, we have $$ -(\partial_i X^i)(p)=(\partial_i\delta_p)(X^i), $$ $$ \Gamma^{j}_{ji}(p)X^i(p)=\Gamma^{j}_{ji}\delta_p(X^i), $$ from which we get $$ (\nabla\delta_p)_i\equiv\nabla_i\delta_p=\partial_i\delta_p-\Gamma^{j}_{ji}\delta_p, $$ which is indeed the rule for covariantly differentiating densities.
