I will abuse notation a little, but I hope you don't find it terrible. After all, I only abuse notation: not children!
Trick 1: Parametrize by Proper Length. We will pick for our affine parameter $\lambda=s$ the proper length. Then the stress energy tensor becomes
$$\tag{1}T^{\alpha\beta}(x)=m\int_{\gamma}u^{\alpha} u^{\beta}\frac{\delta^{(4)}\bigl(x,z(s)\bigr)}{\sqrt{|g|}}\,\mathrm{d}s$$
where $u^{\alpha}=\mathrm{d}x^{\alpha}/\mathrm{d}s$ and $g=\det{g_{\mu\nu}}$.
Trick 2: Covariant Derivative Trick. We can write
$$\nabla_{\mu}f^{\mu}=\frac{1}{\sqrt{|g|}}\partial_{\mu}(\sqrt{|g|}f^{\mu})$$
for arbitrary $f^{\mu}$.
Exercise: Using Poisson's notation (13.2), we have
$$\delta(x,x') = \frac{\delta^{(4)}(x-x')}{\sqrt{|g|}} = \frac{\delta^{(4)}(x-x')}{\sqrt{|g'|}}$$
and thus using our Covariant Derivative trick, find
$$\nabla_{\mu}\delta(x,x')=???$$
This will tell you that
$$\int u^{\alpha}u^{\beta}\nabla_{\alpha}\delta(x,x')\,\mathrm{d}s = \mbox{boundary terms}$$
and thus we can ignore it.
Remark 1. You are in error writing
\begin{alignedat}{1}\nabla_{\beta}T^{\alpha\beta} & =m{\displaystyle \int_{\gamma}\nabla_{\beta}\left[\frac{g^{\alpha}{}_{\mu}g^{\beta}{}_{\nu}\dot{z}^{\mu}\dot{z}^{\nu}}{\sqrt{-g_{\alpha\beta}\dot{z}^{\alpha}\dot{z}^{\nu}}}\delta_{4}\left(x,z\right)\right]d\lambda}=\\
& =m{\displaystyle \int_{\gamma}\nabla_{\beta}\left[\frac{g^{\alpha}{}_{\mu}g^{\beta}{}_{\nu}\dot{z}^{\mu}\dot{z}^{\nu}}{\sqrt{-g_{\alpha\beta}\dot{z}^{\alpha}\dot{z}^{\nu}}}\right]\delta_{4}\left(x,z\right)d\lambda-m{\displaystyle \int_{\gamma}\frac{g^{\alpha}{}_{\mu}g^{\beta}{}_{\nu}\dot{z}^{\mu}\dot{z}^{\nu}}{\sqrt{-g_{\alpha\beta}\dot{z}^{\alpha}\dot{z}^{\nu}}}\nabla_{\nu}\left[\delta_{4}\left(x,z\right)\right]d\lambda,}}
\end{alignedat}
This should have been a simple application of the product rule. That is, the minus sign should be a plus sign.
Remark 2. Why should we expect the right hand side of $\nabla_{\beta}T^{\alpha\beta}=0$? Well, because using Einstein's field equation it's $\nabla_{\beta}G^{\alpha\beta}$ and this is identically zero.
This is why we set
$$\tag{2}m\int_{\gamma}\nabla_{\beta}\left[\frac{g^{\alpha}{}_{\mu}g^{\beta}{}_{\nu}\dot{z}^{\mu}\dot{z}^{\nu}}{\sqrt{-g_{\alpha\beta}\dot{z}^{\alpha}\dot{z}^{\nu}}}\right]\delta\bigl(x,z\bigr)\,\mathrm{d}\lambda=0.$$
...which is precisely the geodesic equation for a point particle as discussed in Poisson's article section 3.
Edit
We can rewrite (2) since ${g^{\alpha}}_{\beta}={\delta^{\alpha}}_{\beta}$ is the Kronecker delta. So
$$\tag{3}m\int_{\gamma}\nabla_{\nu}\left[\frac{\dot{z}^{\alpha}\dot{z}^{\nu}}{\sqrt{-g_{\alpha\beta}\dot{z}^{\alpha}\dot{z}^{\beta}}}\right]\delta\bigl(x,z\bigr)\,\mathrm{d}\lambda=0.$$
But if we pick the arclength as the parameter, this becomes simply
$$\tag{4}m\int_{\gamma}\nabla_{\nu}(u^{\alpha}u^{\nu})\delta\bigl(x,z\bigr)\,\mathrm{d}\lambda=0.$$
Great, but really is
$$\tag{5}\nabla_{\nu}(u^{\alpha}u^{\nu})=0?$$
Lets recall for a geodesic using arclength parametrization we have
$$u_{\mu}u^{\mu}=1\implies u_{\mu}\nabla_{\nu}u^{\mu}=0.$$
Thus (5), when contracted by a non-negative vector (say $u_{\alpha}$) becomes
\begin{alignedat}{1}u_{\alpha}\nabla_{\nu}(u^{\alpha}u^{\nu}) &= u_{\alpha}\underbrace{u^{\nu}\nabla_{\nu}u^{\alpha}}_{=0} + \underbrace{u_{\alpha}u^{\alpha}}_{=1}\nabla_{\nu}u^{\nu}\\
&=\nabla_{\nu}u^{\nu}\end{alignedat}
But this is a continuity-type equation (and if you use trick 2, it really resembles electromagnetism's continuity equation!).
Now we can go back, and by inspection we find
$$\nabla_{\nu}(u^{\alpha}u^{\nu})=\left(\begin{array}{c}\mbox{Geodesic}\\
\mbox{Equation}\end{array}\right)+\left(\begin{array}{c}\mbox{Continuity}\\
\mbox{Equation}\end{array}\right).$$
This is, of course, $\nabla_{\mu}T^{\mu\nu}$. Why should we expect it to be zero?
Well, if the Einstein field equations hold, then
$$G^{\mu\nu}-\kappa T^{\mu\nu}=0$$
and moreover
$$\nabla_{\mu}(G^{\mu\nu}-\kappa T^{\mu\nu})=0.$$
However, $\nabla_{\mu}G^{\mu\nu}=0$ identically thanks to geometry.
Best Answer
It is easier to make this calculation covariantly.
Given a time-like vector field $n^\mu$ representing the direction normal to the hyper-surfaces (in the Friedmann case the homogeneous and isotropic hyper-surfaces), you can always decompose the energy momentum tensor $T_{\mu\nu}$ as $$T_{\mu\nu} = \rho n_\mu n_\nu + 2n_{(\mu}q_{\nu)} + p \gamma_{\mu\nu}+ \Pi_{\mu\nu},$$ where $\rho \equiv T_{\mu\nu}n^\mu n^\nu$ is the energy density, $q_\alpha \equiv - T_{\mu\nu} n^\mu \gamma^\nu{}_\alpha$ the energy (heat) flow, $p \equiv T_{\mu\nu}\gamma^{\mu\nu} / 3$ the pressure and $\Pi_{\mu\nu} \equiv T_{\alpha\beta}\gamma^{\alpha}{}_\mu\gamma^{\beta}{}_{\nu} - p\gamma_{\mu\nu}$ all measured in the hyper-surfaces described by $n^\mu$, $\gamma_{\mu\nu} \equiv g_{\mu\nu} + n_\mu n_\nu$ is the spatial projector/metric and the signature is $(-1,1,1,1)$, i.e., $n_\mu n^\mu = -1$.
The energy momentum tensor "conservation" gives $\nabla^\mu T_{\mu\nu} = 0$. To project this result in the time direction (that is want you want when you put $\nu = 0$) you just multiply this expression by $n^\nu$, i.e., $$n^\nu\nabla^\mu T_{\mu\nu} = \nabla^\mu (T_{\mu\nu}n^\nu) - T_{\mu\nu}\nabla^\mu n^\nu = \nabla^\mu(-\rho n_\mu - q_\mu) - \mathcal{K}^{\mu\nu}T_{\mu\nu},$$ where $\mathcal{K}_{\mu\nu} \equiv \nabla_\mu n_\nu$ is the extrinsic curvature.
Now for a Friedmann metric the extrinsic curvature is simply $\mathcal{K}_{\mu\nu} = \Theta \gamma_{\mu\nu} / 3$ where $\Theta$ is the expansion factor and the heat flow $q_\mu$ is zero, hence, we obtain $$n^\nu\nabla^\mu T_{\mu\nu} = -\left[\dot{\rho}+\Theta(\rho + p) + \nabla_\mu q^\mu + \sigma^{\mu\nu}\Pi_{\mu\nu}\right],$$ where $\dot{\rho} \equiv n^\mu\nabla_\mu \rho$ and $\sigma^{\mu\nu}$ is the traceless symmetric part of the extrinsic curvature which is called shear tensor. At this point one can also introduces a scale factor satisfying $3\dot{a}/a = \Theta$.
The key ingredient to make this calculation correctly is the explicit introduction of the time vector $n^\mu$, note that it is its commutation with the covariant derivative and its derivative in terms of the expansion factor that introduces the missing terms in your result. If you want to make this calculation in the old fashion coordinate way you can just calculate everything in a specific coordinate system with all partial derivatives, Christoffel symbols, etc, and then project in the chosen direction.