Since the space of scalar densities/top forms is one-dimensional, we can make things easier by just computing the component $\nabla_\nu \rho_{1\ldots n}$ instead of having all those extra indices flying around everywhere. Using your formula, this becomes
$$ \nabla_\nu \rho_{1 \ldots n} = \partial_\nu \rho - \rho \left( C^\sigma_{\nu 1} \pi_{\sigma 2 \ldots n} + \cdots + C^\sigma_{\nu n} \pi_{1\ldots (n-1)\sigma}\right).$$
Note that $\pi_{\sigma 2 \ldots n} = \delta_{\sigma}^{1}$, so the first term in the parentheses is $C^\sigma_{\nu 1} \delta_{\sigma}^{1} = C^1_{\nu 1}$. Similarly, if we replace some index $k \in \{ 1, \ldots, n\}$ in $\pi_{1\ldots n}$ with $\sigma$ and contract it with $C^\sigma_{\nu k}$, we get $C^k_{\nu k}$ (no sum); and thus adding these all up we get $$\nabla_\nu \rho_{1 \ldots n} = \partial_\nu\rho -
C^\mu_{\nu \mu} \rho$$ as desired.
If you know how to commute covariant derivatives of an arbitrary tensor, then this is not so hard: since the metric commutes with covariant derivatives, we have $$\begin{align}
\nabla_k \Delta T - \Delta \nabla_k T &= g^{ij}(\nabla_k \nabla_i \nabla_j T - \nabla_i \nabla_j \nabla_k T) \\
&= g^{ij}([\nabla_k, \nabla_i] \nabla_j T +\nabla_i\nabla_k\nabla_j T- \nabla_i \nabla_j \nabla_k T) \\
&= g^{ij}([\nabla_k, \nabla_i] \nabla_j T + \nabla_i([\nabla_k,\nabla_j] T))\tag{1}.\end{align}$$
Here I'm using the abstract index notation, so e.g. $\nabla_i \nabla_j \nabla_k T = (\nabla^3 T)(\partial_i, \partial_j, \partial_k, \ldots)$ is the covariant third derivative.
Since commutators of covariant derivatives always end up just being some contractions with the curvature, we can immediately see from $(1)$ that the formula will have the form $$[\nabla,\Delta]T = R* \nabla T + \nabla( R * T) = R * \nabla T + \nabla R * T \tag{2},$$ where the notation $S * T$ denotes an arbitrary linear combination of contractions of $S \otimes T$. Thus you should expect not only curvature terms, but also derivatives of the curvature. (That is, unless $n = 0$ - meaning $T$ is just a scalar - in which case the second term vanishes because the Hessian of a scalar is symmetric.) For many purposes
$(2)$ is all you need.
To get an explicit expansion of $(1)$ in terms of the curvature tensor, we first need to write down an expression for $[\nabla_k, \nabla_j] T$.$\def\t{\otimes}$ When we take the second covariant derivative of $T = T_{i_1 \ldots i_n}dx^{i_1}\t\cdots\t dx^{i_n}$ and fully expand using the product rule, all the cross terms will vanish thanks to the antisymmetrization; so the derivative commutator also satisfies a Leibniz rule: we find $$[\nabla_k, \nabla_j] T = [\nabla_k,\nabla_j](T_{i_1 \ldots i_n})dx^{i_1}\t\cdots\t dx^{i_n}\\ + T_{i_1 \ldots i_n} ([\nabla_k,\nabla_j]dx^{i_1})\t\cdots \t dx^{i_n} \\ \vdots \\+ T_{i_1 \ldots i_n} dx^{i_1}\t\cdots \t([\nabla_k,\nabla_j]dx^{i_n}).$$
Since the curvature tensor is characterized by $R_{kj}{}^i{}_l dx^l = [\nabla_k,\nabla_j]dx^i$ and the first two derivatives of a scalar commute, this becomes $$[\nabla_k, \nabla_j] T=T_{i_1 \ldots i_n}(R_{kj}{}^{i_1}{}_ldx^l \t \cdots \t dx^{i_n} + \cdots + R_{kj}{}^{i_n}{}_l dx^{i_1} \t \cdots \t dx^l),$$ which (doing some index gymnastics) we can write as
$$[\nabla_k,\nabla_j]T_{i_1\ldots i_n} = R_{kj}{}^l{}_{i_1} T_{l i_2 \ldots i_n}+\cdots + R_{kj}{}^l{}_{i_n}T_{i_1 \ldots i_{n-1} l}.$$
Since you can now commute derivatives of an arbitrary $n$-tensor, you can do it for an $n+1$-tensor too; so if you really want to write out the full mess you can plug this in to $(1)$ for the two commutator terms and then expand everything with the product rule. I don't think there's much insight to be gained by typing it all out - it's just going to be a mess of contractions, some of them of $R$ with $\nabla T$ and some of $\nabla R$ with $T$. Perhaps try it yourself for a very small value of $n$.
Best Answer
I did the case for a $(0, 1)$ tensor $A = A_i$, it's
$$ (\Delta A)_i = \Delta (A_i) - g^{jk}\partial _k \Gamma_{ij}^l A_l - 2 g^{jk}\Gamma_{ij}^l A_{l,k} +g^{jk}\Gamma_{ik}^l \Gamma_{jl}^m A_m + g^{jk} \Gamma_{jk}^l \Gamma_{il}^m A_m$$
In general there is no hope that $(\Delta A)_i = \Delta (A_i)$, since the equations should be all coupled and the term $g^{jk}\partial _k \Gamma_{ij}^l A_l$ involves the curvature. If your manifold is something explicit, might be you can simplify the term though.