I want to show
$$
(\delta S)Z = \delta (S(\cdot, Z))+\frac{1}{2}\langle S, \mathcal L_Zg\rangle
\tag{1}
$$
where $\delta= -tr_{12}\nabla$ is divergence operator, $\mathcal L$ is Lie derivative, $S\in\Gamma(Sym^2 T^*M)$
What I try: I use abstract index notation. First $$ (\delta S)Z =-(\nabla^a S_{aj})Z = -(\nabla^a S_{aj})Z^j $$ Besides, $$ \delta (S(\cdot,Z)) =-\nabla^a[S(\cdot, Z)]_a $$ since $S(\cdot, Z)= S_{ij}Z^j$, I think $[S(\cdot, Z)]_a=S_{aj}Z^j$. Therefore $$ \delta(S(\cdot, Z)) =-\nabla^a(S_{aj}Z^j)=-(\nabla ^a S_{aj})Z^j - S_{aj}(\nabla^a Z^j) =(\delta S)Z- S_{aj}(\nabla^a Z^j) \tag{2} $$ On the other hand, when $\delta$ is restricted on $\Gamma (Sym^2 T^*M)$, its formal adjoint is $\omega\rightarrow \frac{1}{2}\mathcal L_{\omega^\sharp}g$. Therefore, I have $$ \frac{1}{2}\langle S, \mathcal L_Z g \rangle = \langle \delta(S), Z^\flat \rangle = \langle -\nabla^a S_{aj} , Z^i {g_{ik}} \rangle = -(\nabla^a S_{aj}) Z^j \tag{3} $$ However, I can't get (1) from (2) and (3).
If the answer is complex, a picture of handwritten draft is enough for me. Thanks very much.
It's good that you are trying to use the abstract index notation, but you need one more piece for you puzzle, the formula for the Lie derivative of the metric: $$ (\mathcal L_Z g)_{a b} = \nabla_a Z_b + \nabla_b Z_a = 2 \nabla_{(a} Z_{b)} $$ where $Z^a$ is a vector field, and $Z_a = g_{a b} Z^b$ is the corresponding $1$-form. Of course, $\nabla$ is the Levi-Civita connection, corresponding to the Riemannian metric $g_{a b}$, and we use the latter along with its inverse $g^{a b}$ to raise and lower the indices without mention (the Ricci calculus).
With that in mind, the calculation will be straightforward: $$ \delta (S(\cdot, Z)) = - \nabla^a (S_{a b} Z^b) = - (\nabla^a S_{a b})Z^b - S_{a b} \nabla^a Z^b = \\ (\delta S) Z - S_{a b} \nabla^{(a} Z^{b)} = (\delta S) Z - \tfrac{1}{2} S_{a b} (\mathcal L_Z g)^{a b} = (\delta S) Z - \tfrac{1}{2} \langle S, \mathcal L_Z g \rangle $$
Here I use the notation for the symmetric part $t_{(a b)} := \tfrac{1}{2}(t_{a b} + t_{b a}) $ of a $2$-tensor $t_{a b}$ and the fact that if a tensor $s_{a b}$ is symmetric ($s_{a b} = s_{b a} = s_{(a b)} $), then for any tensor $t_{a b}$ we have the identity: $$ s_{a b} t^{a b} = s_{(a b)} t^{a b} = s_{a b} t^{(a b)} $$