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.
The geodesic equation can be written equivalently when $X$ is parametrized with proper time:
$$
\frac{d^2X^\mu}{d\tau^2} + \Gamma^\mu_{\nu\alpha}\frac{dX^\nu}{d\tau} \frac{dX^\alpha}{d\tau} = 0\,.
$$
Here, $\tau\mapsto X^\mu(\tau)$ is a world line whose four velocity is $\frac{dX^\mu}{d\tau}$ and whose acceleration is $\frac{d^2X^\mu}{d\tau^2}\,.$
Now suppose we have a four-velocity field $U^\mu$ at every point of space time (or at least around the world line $X^\mu$) such that
$$
\frac{dX^\mu(\tau)}{d\tau}=U^\mu(X^\mu(\tau))\,.
$$
By the chain rule
$$\tag{1}
\frac{d^2X^\mu(\tau)}{d\tau^2}=\frac{dX^\mu(\tau)}{d\tau}\partial_\nu U^\mu(X^\mu(\tau))=(U^\nu\partial_\nu U^\mu)(X^\mu(\tau))\,.
$$
Dropping the clumsy argument $X^\mu(\tau)$ the geodesic equation then implies
$$
U^\nu\partial_\nu U^\mu+\Gamma^\mu_{\nu\alpha}U^\nu U^\alpha=U^\nu\nabla_\nu U^\mu=0.
$$
Four momentum $P^\mu$ is just rest mass $m_0$ times four velocity $U^\mu\,.$ This does not alter the geodesic equation.
Few remarks.
In curved space time $\frac{d^2X^\mu}{d\tau^2}$ is not the four acceleration. Since by (1) it is $U^\nu\partial_\nu U^\mu$ which does not transform as a tensor due to the partial derivative there.
The correct way of defining four acceleration is
$$\tag{2}
A^\mu\equiv\frac{d^2X^\mu}{d\tau^2}+\Gamma^\mu_{\nu\alpha}\frac{dX^\nu}{d\tau}\frac{dX^\alpha}{d\tau}\,.
$$
As above this is seen to be equal to $U^\nu\nabla_\nu U^\mu\,.$
This is the covariant derivative of $U$ in the direction of $U\,.$
Some authors also write $\dot U^\mu$ for $A^\mu$
and the shortest form of the geodesic equation is
$$
\boxed{A^\mu=0\,.}
$$
Best Answer
Comments to the question (v1):
Note first of all that there isn't a canonical/unique prescription to assign values to $\tau_i$ and $\tau_f$ for a given path $\gamma$.
In particular, it is not assumed that $\tau_i$ and $\tau_f$ are kept fixed during the variation.
In contrast, the parameter endpoints $\lambda_i$ and $\lambda_f$, and the path endpoints $\gamma_i$ and $\gamma_f$, are kept fixed during the variation.
Only the difference $\Delta\tau = \tau_f - \tau_i$ is important. In fact $$\tag{1}\Delta\tau~=~\int_{\lambda_i}^{\lambda_f} \!d\lambda ~\sqrt{-g_{\mu\nu}(\gamma(\lambda)) ~\dot{\gamma}^{\mu}(\lambda)~\dot{\gamma}^{\nu}(\lambda)}$$ is the proper time of the path and precisely the functional that we want to vary.
This functional (1) does depend on the path $\gamma$.