Why are the Klein-Gordon equations warranted from the conservation of the energy-momentum tensor?

If we have an action with a scalar field non-minimally coupled to the gravity: $$\int dx^4 \sqrt{-g}(-\frac{1}{2}\partial_\mu\phi\partial^\mu\phi-\frac{1}{2}\zeta R\phi^2-V(\phi)).....(1)$$ varying respect to the scalar field, we can get a Klein-Gordon like equation: $$\nabla^\mu\nabla_\mu \phi=\zeta R \phi+V'(\phi)....(2)$$ And they say "The equation (2) is warranted from the conservation of the energy-momentum tensor of the scalar field due to the diffeomorphism invariance of action. $$\nabla^\mu T_{\mu\nu}=0$$

Why is this happening?

Paper where I see this

The demonstration compounds of two pieces: first you define the energy-momentum tensor of the matter fields in the nonvacuum Einstein's field equations, secondly apply the diffeomorphism invariance to the complete action.

Define the energy-momentum tensor of the matter fields
$$S = \frac{1}{16 \pi} S_H + S_M$$
where:
$$c = G = 1$$ natural units
$$S$$ complete action for gravity coupled to a set of matter fields
$$S_H$$ Hilbert action
$$S_H = \int \sqrt{-g} R d^n x$$
$$S_M$$ action for matter
Applying the principle of least action by varying the complete action with respect to the metric, we get:
$$\frac{1}{\sqrt{-g}} \frac{\delta S}{\delta g^{\mu \nu}} = \frac{1}{16 \pi} (R_{\mu \nu} - \frac{1}{2} R g_{\mu \nu}) + \frac{1}{\sqrt{-g}} \frac{\delta S_M}{\delta g^{\mu \nu}} = 0$$
$$T_{\mu \nu} = -2 \frac{1}{\sqrt{-g}} \frac{\delta S_M}{\delta g^{\mu \nu}}$$ energy-momentum tensor (definition)
Hence:
$$R_{\mu \nu} - \frac{1}{2} R g_{\mu \nu} = 8 \pi T_{\mu \nu}$$ Einstein's field equations

Apply the diffeomorphism invariance to the complete action
$$S = \frac{1}{16 \pi} S_H [g_{\mu \nu}] + S_M [g_{\mu \nu}, \psi^i]$$
where:
$$g_{\mu \nu}$$ metric tensor
$$\psi^i$$ matter fields
The Hilbert action $$S_H$$ is diffeomorphism invariant, that is coordinate invariant, when considered in isolation. As the complete action is to be diffeomorphism invariant as well, the matter action $$S_M$$ must also be.
Under a diffeomorphism:
$$\delta S_M = \int d^n x \frac{\delta S_M}{\delta g_{\mu \nu}} \delta g_{\mu \nu} + \int d^n x \frac{\delta S_M}{\delta \psi^i} \delta \psi^i$$
Since the gravitational part of the complete action does not involve the matter fields, the variation of $$S_M$$ with respect to $$\psi^i$$ will vanish. Therefore also the variation of $$S_M$$ with respect to $$g_{\mu \nu}$$ must vanish.
If the diffeomorphism is generated by an infinitesimal vector field $$V^\mu (x)$$, the infinitesimal change in the metric is given by the Lie derivative along $$V^\mu$$.
$$\delta g_{\mu \nu} = \mathcal L_V g_{\mu \nu} = 2 \nabla _{( \mu} V_{\nu )}$$
Setting $$\delta S_M = 0$$ implies:
$$0 = \int d^n x \frac{\delta S_M}{\delta g_{\mu \nu}} \nabla_\mu V_\nu = - \int d^n x \sqrt{-g} V_\nu \nabla_\mu (\frac{1}{\sqrt{-g}} \frac{\delta S_M}{\delta g_{\mu \nu}})$$
The symmetrization of $$\nabla _{( \mu} V_{\nu )}$$ was dropped as $$\frac{\delta S_M}{\delta g_{\mu \nu}}$$ is already symmetric.

Law of energy-momentum conservation
Demanding that the last equation holds for diffeomorphisms generated by arbitrary vector fields $$V^\mu$$ and using the definition of the energy-momentum tensor, we get the law of energy-momentum conservation:
$$\nabla_\mu T^{\mu \nu} = 0$$
Note: The derivation assumes that no matter fields appear in the gravitational action.