Show that $(divR) (X, Y,Z) = (\nabla_X Ric) (Y,Z) − (\nabla_Y Ric) (X,Z).$

In course of solving Riemannian Geometry By Peter Petersen Chap. 2, I stuck on the following problem:

Show that in a Riemmanian manifold if $R$ is the $(1, 3)$ curvature tensor and $Ric$ the $(0, 2)$ Ricci tensor, then $(div~R) (X, Y,Z) = (\nabla_X Ric) (Y,Z) − (\nabla_Y Ric) (X,Z).$

I am absolutely clueless. Can anyone help me to solve it?

Consider an orthonormal frame $\{E_i\}.$ Consider $X,Y,Z$ such that at a given point $p$ it is $\nabla_{E_i}E_j=\nabla_{E_i}X=\nabla_{E_i}Y=\nabla_{E_i}Z=0.$ Then, at $p,$ by definition,
$$(\mathrm{div} R)(X,Y,Z) = \sum_{i=1}^n g((\nabla_{E_i}R)(X,Y,Z),E_i).$$ Now, having in mind the second Bianchi identity, $$(\nabla_UR)(V,W)+(\nabla_VR)(W,U)+(\nabla_WR)(U,V)=0,$$ one has
$$g((\nabla_{E_i}R)(X,Y,Z),E_i)=-g((\nabla_{X}R)(Y,E_i,Z),E_i)-g((\nabla_{Y}R)(E_i,X,Z),E_i).$$ But, since $\nabla$ commutes with traces, it is
$$\sum_{i=1}^ng((\nabla_{Y}R)(E_i,X,Z),E_i)=(\nabla_YRc)(X,Z)$$ and $$\sum_{i=1}^ng((\nabla_{X}R)(E_i,Y,Z),E_i)=-(\nabla_XRc)(Y,Z).$$ So, we have obtained the desired equality.