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 gravity,
$$\int dx^4 \sqrt{-g}(-\frac{1}{2}\partial_\mu\phi\partial^\mu\phi-\frac{1}{2}\zeta R\phi^2-V(\phi)), \tag{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) \tag{2}.$$
arXiv: 1506.03545 [gr-qc] says

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?
 A: 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.
