What is the purpose of emphasizing that an action is invariant under diffeomorphism?
After checking the book 'Differential Geometry and Lie Groups for Physicists' by Marian Fecko in section 16.4.1, I believe that I am close to understanding what physicists mean by invariance of an action under diffeomorphism. In what follows, I explain Marian Fecko's discussion on diffeomorphism invariance of an action.
Consider a classical field theory $\phi$ on a Riemannian manifold $(M,g)$, where $g_{ab}$ is its metric. We say the action natural with respect to diffeomorphism in the following sense.
Let us define such a differential form $\Omega[\phi,g]$ defined via the action $$S[\phi,g]=\int_{D}\Omega[\phi,g]=\int_{D}L(\phi,\nabla\phi,g)\omega_{g}$$ where $D\subset M$ is a submanifold and $\omega_{g}$ is the volume form on $D$ associated with the metric $g$.
Let us then assume that the variation of fields and the variation of spacetime coordinates vanish on the boundary $\partial D$. In mathematical terminology, this requirement corresponds to a flow $\Psi_{t}:M\rightarrow M$ that is arbitrary inside $D$ but vanishes on $\partial D$.
We say that the action is invariant (or natural) under the diffeomorphism $\Psi_{t}$ if the pull-back satisfies $$\Psi_{t}^{\ast}(\Omega[\phi,g])=\Omega[\Psi_{t}^{\ast}(\phi),\Psi_{t}^{\ast}(g)].$$
Under such a requirement, since the flow does not move points on the boundary, under any infinitesimal amount of variation $\Psi_{\epsilon}$ generated by a vector field $V$ in $M$, we have $\Psi_{\epsilon}(D)=D$. It follows that $$S[\phi,g]=\int_{D}\Omega[\phi,g]=\int_{\Psi_{\epsilon}(D)}\Omega[\phi,g]=\int_{D}\Psi_{\epsilon}^{\ast}(\Omega[\phi,g])$$ $$=\int_{\Psi_{\epsilon}(D)}\Psi_{\epsilon}^{\ast}(\Omega[\phi,g])=\Psi_{\epsilon}^{\ast}(S[\phi,g])$$ $$=\int_{D}\Omega[\Psi_{\epsilon}^{\ast}(\phi),\Psi_{\epsilon}^{\ast}(g)]=\int_{D}\Omega[\phi+\epsilon\mathcal{L}_{V}\phi,g+\epsilon\mathcal{L}_{V}g]+o(\epsilon),$$ where $\mathcal{L}_{V}$ is the Lie derivative along the flow, and we have used integration by substitution in the first line. Since classical fields must be on-shell, $\phi$ extremizes the action $S[\phi,g]$. We then have $$\int_{D}\Omega[\Psi_{\epsilon}^{\ast}(\phi),\Psi_{\epsilon}^{\ast}(g)]=\int_{D}\Omega[\phi+\epsilon\mathcal{L}_{V}\phi,g+\epsilon\mathcal{L}_{V}g]=\int_{D}\Omega[\phi,g+\epsilon\mathcal{L}_{V}g].$$ and so $$S[\phi,g]=\int_{D}\Omega[\phi,g]=\Psi_{\epsilon}^{\ast}(S[\phi,g])=\int_{D}\Omega[\phi,g+\epsilon\mathcal{L}_{V}g]+o(\epsilon).$$ Therefore, by definition of energy-momentum tensor, we have $$S[\phi,g]=S[\phi,g]-\epsilon\int_{D}\frac{1}{2}(\mathcal{L}_{V}g)_{ab}T^{ab}\omega_{g}+o(\epsilon),$$ and so $$\int_{D}(\mathcal{L}_{V}g)_{ab}T^{ab}\omega_{g}=0,$$ for arbitrary variation $\delta g$, we conclude that the energy-momentum is conserved i.e. $\nabla_{a}T^{ab}=0$.
In summary, the diffeomorphism invariance of an action is a crucial condition for the conservation of energy-momentum tensor of the system.
I hope this can contribute the understanding of those who once were also confused by it. Welcome to clarify any mistakes and any misunderstandings I have.
Perhaps a simple example is in order: The action for a non-relativistic free particle
$$ S[x]~=~\int \! dt ~L, \qquad L ~=~ \frac{m}{2}\dot{x}^2, $$
is not form invariant under time-reparametrizations
$$t\quad\longrightarrow \quad t^{\prime} ~=~f(t). $$
In contrast, modern fundamental physics (such as, e.g. string theory) is believed to be geometric, and the action formulation are expected to be reparametrization and diffeomorphism invariant.