The Lie derivative is defined as picture below.$X,Y,Z$ is vector fields,and $g$ is Riemannian metric. I try to compute $$ \mathcal L_Xg(Y,Z)=X(g(Y,Z))-g(\mathcal L_XY,Z)-g(Y,\mathcal L_XZ) $$
Firstly, \begin{align} \mathcal L_Xg(Y,Z) &=(\frac{d}{dt}(*\varphi_t)g(Y,Z))|_{t=0}\\ &=\frac{d}{dt}|_{t=0}g_{\varphi_t}(\varphi_t^*(Y),\varphi^*_t(Z)) \end{align}
Then, I don't know how to do it.
To get the identity you want, you don't need the definition of the Lie derivative, the properties listed in Lemma 1.6 of your post are sufficient. The main issue here is that you can consider $g\otimes Y\otimes Z$, which is a section of $\otimes ^2T^*M\otimes\otimes^2 TM$, whose complete contraction equals $g(Y,Z)\in C^\infty(M,\mathbb R)$. By 1. of Lemma 1.6, the Lie derivative of this complete contraction in the direction of $X$ equals $X\cdot g(Y,Z)$. Since the complete contraction is a composition of traces, you can compute this as the complete contraction of $$\mathcal L_X(g\otimes Y\otimes Z)=(\mathcal L_Xg)\otimes Y\otimes Z+g\otimes[X,Y]\otimes Z+g\otimes Y\otimes[X,Z]. $$ (Here part 3. of Lemma 1.6 was used to expand the Lie derivative and part 2. was applies to compute $\mathcal L_X Y$ and $\mathcal L_X Z$.) Now the complete contraction of the right hand side is just $(\mathcal L_X g)(Y,Z)+g([X,Y],Z)+g(Y,[X,Z])$ and bringing the last two terms to the other side, the result follows.