Consider here a background metric $g_{\mu\nu}$, we impose a perturbation $g_{\mu\nu}+\epsilon h_{\mu\nu}$ with $\epsilon\ll1$. Then we can write down the modified Einstein-Hilbert action with zero cosmological constant: \begin{align} S&=\int\mathrm{d}^{D}x\sqrt{-g}\left(1-\epsilon\frac{1}{2}h_{\mu\nu}g^{\mu\nu}\right)\left(R^{0}+\epsilon R^{1}\right)\\ &=S^{0}+\epsilon\int\mathrm{d}^{D}x\sqrt{-g}\left(R^{1}-\frac{1}{2}R^{0}g_{\mu\nu}h^{\mu\nu}\right)+\mathcal{O}(\epsilon^2)\, \end{align} where $R^{0}$ is the Ricci scalar w.r.t $g_{\mu\nu}$, and $R^{1}$ is the correction of $R$ at the leading order which can be found in this answer.
To my knowledge, if for any: $$\delta_{\zeta}h_{\mu\nu}=\nabla_{\mu}\zeta_{\nu}+\nabla_{\nu}\zeta_{\mu}=\mathcal{L}_{\zeta}g_{\mu\nu} \, \ \ \ \delta_{\zeta}g_{\mu\nu}=0 \,$$ the action is invariant at each order, this theory is gauge invariant, and the perturbation field $h_{\mu\nu}$ is a spin-2 gauge field. For example, when the background is flat $g_{\mu\nu}=\eta_{\mu\nu}$, $R^{0}=0$, one can check this easily, for an arbitrary $\zeta \ $: \begin{align} \delta_{\zeta}S&=\delta_{\zeta}S^0+\epsilon\int\mathrm{d}^{D}x\sqrt{-g}\delta_{\zeta}R^{1} +\mathcal{O}(\epsilon^2)\\ &=0. \end{align} But for a curved background $R^{0}\ne 0$, in principle one should show whether: $$\int\mathrm{d}^Dx\sqrt{-g}\left(\delta_{\zeta}R^{1}-\frac{1}{2}R^{0}g_{\mu\nu}\mathcal{L}_{\zeta}g^{\mu\nu}\right)=0 \ .$$ According to the P4 eq(2.4) in this paper 2110.13512, one find \begin{equation} \delta_{\zeta}R^{1}=\mathcal{L}_{\zeta}R^{0}. \end{equation} Thus, for any $g_{\mu\nu}$ \begin{align} \mathcal{L}_{\zeta}R-\frac{1}{2}Rg_{\mu\nu}\mathcal{L}_{\zeta}g^{\mu\nu}&=\mathcal{L}_{\zeta}(g^{\mu\nu}R_{\mu\nu})-\frac{1}{2}Rg_{\mu\nu}\mathcal{L}_{\zeta}g^{\mu\nu}\\ &=g^{\mu\nu}\mathcal{L}_{\zeta}R_{\mu\nu}+\left[R_{\mu\nu}-\frac{1}{2}Rg_{\mu\nu}\right]\mathcal{L}_{\zeta}g^{\mu\nu}\\ &=g^{\mu\nu}\mathcal{L}_{\zeta}R_{\mu\nu}+\mathrm{EoM}\mathcal{L}_{\zeta}g^{\mu\nu} \end{align} if $g_{\mu\nu}$ is a solution of field equation, the gauge invariance of a linearized theory if and only if $\delta_{\zeta}R^1_{\mu\nu}=\mathcal{L}_{\zeta}R^0_{\mu\nu}=0$, as expected. But it is not easy to calculate. So on the other hand, we express it as \begin{align} \mathcal{L}_{\zeta}R-\frac{1}{2}Rg_{\mu\nu}\mathcal{L}_{\zeta}g^{\mu\nu}&=\mathcal{L}_{\zeta}R+\frac{1}{2}Rg^{\mu\nu}\mathcal{L}_{\zeta}g_{\mu\nu}\\ &=\zeta^\mu\nabla_\mu R+R\nabla_\mu\zeta^\mu\\ &=\nabla_\mu(\zeta^\mu R) \end{align} This equation states that a theory is gauge invariant if and only if $\zeta^\mu R^{0}\big|_{\partial\mathcal{M}}=0$. This result is suspicious to me, is something wrong here?
Especially, for $dS$ spacetime, this boundary term is non-zero. So is there a gauge-invariant linearized theory in $dS$?