Suppose $M$ to be orientable and oriented and let $\mu\in\Omega^m(M)$ be a volume form ($m=\dim M$). Let $\nabla$ be a linear connection on $\tau:TM\rightarrow M$.
Given a vector field $X\in\mathcal D(M)$, the divergence of $X$ with respect to the volume form $\mu$ is $$ \mathrm{Div}_\mu(X):=(\mathscr L_X\mu)/\mu, $$ where $\mathscr L_X$ is the Lie derivative w.r.t. $X$. The divergence of $X$ with respect to the connection $\nabla$ is $$ \mathrm{Div}_\nabla(X):=\mathrm{Tr}(\nabla X). $$
Gauss' theorem always holds for divergence with respect to a volume form. Specifically, if $\Omega\subseteq M$ is a compact $m$ dimensional submanifold with boundary, then $$ \int_\Omega\mathrm{Div}_\mu(X)\mu=\int_\Omega\mathscr L_X\mu=\int_\Omega\mathrm{d}(X\rfloor\mu)=\int_{\partial\Omega}X\rfloor\mu. $$
Now if $(U,x)$ is a local chart for $M$, and if $\mu=\rho\, \mathrm dx^1\wedge\dots\wedge\mathrm dx^m$, then $$ \mathrm{Div}_\mu(X)=\frac{1}{\rho}\frac{\partial}{\partial x^i}\left(\rho X^i\right),\quad X=X^i\frac{\partial}{\partial x^i}. $$
Meanwhile, let us use the convention for covariant differentiation that the index corresponding to the direction of differentiation is the last one, $$ \nabla_k X^i=\frac{\partial X^i}{\partial x^k}+\Gamma^i_{\ jk}X^j. $$ Then $$ \mathrm{Div}_\nabla(X)=\frac{\partial X^i}{\partial x^i}+\Gamma^i_{\ ji}X^j.$$ The covariant derivative of the volume form is $$ \nabla_i\rho=\frac{\partial\rho}{\partial x^i}-\Gamma^j_{\ ji}\rho. $$
We now want to look for a condition where the two notions of divergence coincide. Expanding the divergence with respect to volume form, we get $$ \mathrm{Div}_\mu(X)=\frac{\partial X^i}{\partial x^i}+X^i\frac{\partial}{\partial x^i}\ln\rho=\frac{\partial X^i}{\partial x^i}+\Gamma^j_{\ ji}X^i+\rho^{-1}X^i\nabla_i\rho. $$
Now, the torsion of the connection is $$ T^i_{\ jk}=\Gamma^i_{\ kj}-\Gamma^i_{\ jk}, $$ and taking a trace gives $$ T_j=\Gamma^k_{\ kj}-\Gamma^k_{\ jk},\quad T_j=T^k_{\ jk}, $$ with which we have $$ \Gamma^j_{\ ji}X^i=T_i X^i+\Gamma^j_{\ ij}X^i, $$ which then gives $$ \mathrm{Div}_\mu(X)=\mathrm{Div}_\nabla(X)+\left(T_i+\rho^{-1}\nabla_i\rho\right)X^i. $$
We then obtain that the condition that the two notions of divergence coincide for all vector fields is that the connection and volume form should satisfy the equation $$ T_i=-\rho^{-1}\nabla_i\rho. $$ In this case, a Gauss theorem applies to the divergence $\mathrm{Div}_\nabla$.
When the vector part ($T_i$) of the torsion vanishes (which of course is the case when the torsion altogether vanishes), then this condition reduces to $\nabla_i\rho=0$, or equivalently $\nabla\mu=0$, i.e. the connection must be volume-preserving.