Let $P$ is an arbitrary zero point at which we expand around. Then expanding metric around $P$ is :
$$g_{ij} = \left. g_{ij} \right\vert_{P} + \frac{1}{2!} \left. \left( \partial_k g_{ij} \right) \right\vert_{P} x^k + \frac{1}{3!} \left. \left( \partial_{k_2}\partial_{k_1} g_{ij} \right) \right\vert_{P} x^{k_1} x^{k_2} + \cdots$$
Curvature tensor generally describes $\nabla_{[i}\nabla_{j]}$. Examining the summation of covariant derivatives of curvature tensor, we can find the highest-order derivatives of the metric appearing in the summation above.
From the Gauss Lemma :
$$x_k=g_{kl}x^l$$
Differentiating 1, 2, 3, ... times, we get
$$\delta_{ki} = \partial_i g_{kj} x^j + g_{ji}$$
$$0=\partial_i \partial_l g_{kj} x^j + \partial_i g_{kl} + \partial_l g_{ki} \;\;\; \rightarrow \;\;\; \left. \left( \partial_i g_{kl} + \partial_l g_{ki} \right) \right\vert_{P}=0$$
$$0=\partial_i \partial_l \partial_s g_{kj} x^j + \partial_i \partial_l g_{ks} + \partial_i \partial_s g_{kl} + \partial_l \partial_s g_{ki} \;\;\; \rightarrow \;\;\; \left. \left( \partial_i \partial_l g_{ks} + \partial_i \partial_s g_{kl} + \partial_l \partial_s g_{ki} \right) \right\vert_{P}=0$$
For notational convenience, let $\partial_{i_1} \partial_{i_2} \cdots \partial_{i_n} g_{ij} = (i_1 i_2\cdots i_n \vert ij)$.
We want to make a summation terms like
$$\nabla_{k_n} \cdots \nabla_{k_3} R_{i k_2 j k_1} x^{k_n} \cdots x^{k_1}$$
First look at the curvature tensor :
$$R_{abcd}=\frac{1}{2} \left( \partial_b \partial_c g_{ad} -\partial_a \partial_c g_{bd} - \partial_b \partial_d g_{ac} + \partial_a \partial_d g_{bc}\right) + \text{a polynomial in lower partial derivs of }g$$
for each $i, j, k_1, \cdots, k_n$, its derivatives are (★) :
$$\begin{align*} \left. \nabla_{k_n} \cdots \nabla_{k_3} R_{i k_2 j k_1} \right\vert_{P} & = \frac{1}{2} \left. \left( \partial_{k_n} \cdots \partial_{k_3} \left( \partial_j \partial_{k_2} g_{k_1 i} - \partial_j \partial_i g_{k_1 k_2} - \partial_{k_1} \partial_{k_2} g_{ij} + \partial_{k_1} \partial_i g_{j k_2} \right) \right) \right\vert_{P} \\ & + \text{a polynomial in lower partial derivs of }g \end{align*}$$
(when $n=2, 3$, lower order partial derivs turns out to be zero)
Using our reduced notation, Gauss Leemma states :
$$(i \vert lk) + (l \vert ik) = 0, \;\;\; (il \vert sk) + (ls \vert ik) + (si \vert lk) = 0, \;\;\; \cdots$$
In order to make the pattern for n-th order Gauss Lemma, fix one of these indices and keep adding successive lists permuting $(n+1)$ indices cyclically. Then n-th order Gauss Lemma is expressed as :
$$(i_1 i_2 \cdots i_{n-3} i_{n-2} \vert i_{n-1} i_n) + (i_2 i_3 \cdots i_{n-2} i_{n-1} \vert i_1 i_n) + \cdots + (i_{n-1} i_1 \cdots i_{n-4} i_{n-3} \vert i_{n-2} i_n) = 0$$
For $n=2$, we can find the coefficient of $x^k x^l$ in $\left. R_{ikjl} \right\vert_{P} x^k x^l$ to be :
$$\begin{align*} \left. R_{ikjl} \right\vert_{P} + \left. R_{iljk} \right\vert_{P} & = -\frac{1}{2} \left[ (ji \vert lk) - (jk \vert li) \right] + \frac{1}{2} \left[ (li \vert jk) - (lk \vert ji) \right] -\frac{1}{2} \left[ (ji \vert kl) - (jl \vert ki) \right] + \frac{1}{2} \left[ (ki \vert jl) - (kl \vert ji) \right] \\ & =-\frac{1}{2} \left[ 2(kl \vert ij) - \left\{ (jk \vert li) + (jl \vert ki) \right\} - \left\{ (ik \vert lj) + (il \vert kj) \right\} + 2(ij \vert kl) \right] \end{align*}$$
Likewise, for n=3, the coefficent of $x^{k_1} x^{k_2} x^{k_3}$ is :
$$\nabla_{k_3} R_{i k_1 j k_2} + \nabla_{k_3} R_{i k_2 j k_1} + \nabla_{k_2} R_{i k_1 j k_3} + \nabla_{k_2} R_{i k_3 j k_1} + \nabla_{k_1} R_{i k_2 j k_3} + \nabla_{k_1} R_{i k_3 j k_2}$$
$$= -\frac{1}{2} \left[ 6(k_3 k_2 k_1 \vert ij) - 2\left\{ (k_3 k_2 j \vert k_1 i) + (k_3 k_1 j \vert k_2 i) + (k_2 k_1 j \vert k_3 i) \right\} \right.$$
$$\left. - 2\left\{ (k_3 k_2 i \vert k_1 j) + (k_3 k_1 i \vert k_2 j) + (k_2 k_1 i \vert k_3 j) \right\} + 2\left\{ (k_3 i j \vert k_2 k_1) + (k_2 i j \vert k_3 k_1) + (k_1 i j \vert k_2 k_3) \right\}\right]$$
Then now examining the coefficient of $x^{k_n} \cdots x^{k_1}$ in $\left. \nabla_{k_n} \cdots \nabla_{k_3} R_{i k_2 j k_1} \right\vert_{P}$, for some arrangement of the indices $k_1, \cdots, k_n$, we get (all of distinct $n!$ arrangement summed) :
$$\left. \sum_{\sigma} \nabla_{k_{\sigma(n)}} \cdots \nabla_{k_{\sigma(3)}} R_{i k_{\sigma(2)} j k_{\sigma(1)}} \right\vert_{P}$$
Substituting from (★), this coefficient becomes
$$-\frac{1}{2} \left[ n! (k_n \cdots k_1 \vert i j) - (n-1)! \sum_{l=1}^{n} (k_n \cdots k_{l+1} k_{l-1} \cdots k_1 j \vert k_l i) - (n-1)! \sum_{l=1}^{n} (k_n \cdots k_{l+1} k_{l-1} \cdots k_1 i \vert k_l j) \right.$$
$$\left.+ 2(n-2)! \sum_{l_2 = 1}^{n} \sum_{l_1=1}^{l_2-1} (k_n \cdots k_{l_1+1} k_{l_1-1} \cdots k_{l_2+1} k_{l_2-1} \cdots k_1 j i \vert k_{l_2} k_{l_1}) \right] + \text{a polynomial in lower order partial derivs of }g$$
From Gauss Lemma, we can simplify the 2nd & 3rd terms as :
$$\sum_{l=1}^{n} (k_n \cdots k_{l+1} k_{l-1} \cdots k_1 j \vert k_l i) = -(k_n \cdots k_1 \vert i j)$$
$$\sum_{l=1}^{n} (k_n \cdots k_{l+1} k_{l-1} \cdots k_1 i \vert k_l j) = -(k_n \cdots k_1 \vert i j)$$
Also, sum of all $_{n+2}C_{n}$ tuples like $(a_1 \cdots a_n \vert a_{n+1} a_{n+2})$ is zero. Then we can simplify the 4th term as :
$$\sum_{l_2 = 1}^{n} \sum_{l_1=1}^{l_2-1} (k_n \cdots k_{l_1+1} k_{l_1-1} \cdots k_{l_2+1} k_{l_2-1} \cdots k_1 j i \vert k_{l_2} k_{l_1}) = (k_n \cdots k_1 \vert i j)$$
which is to say that
$$ \begin{align*} \nabla_{k_n} \cdots \nabla_{k_3} R_{i k_2 j k_1} x^{k_n} \cdots x^{k_1} & = -\frac{n!+2(n-1)!+2(n-2)!}{2n!} \left( \partial_{k_n} \cdots \partial_{k_1} g_{ij} \right) x^{k_n} \cdots x^{k_1} \\ & + \left(\text{a polynomial in lower order partial derivs of }g\right) x^{k_n} \cdots x^{k_1} \end{align*}$$
Reference : Niles G. Johnson, 2003, "The Riemannian Metric and Curvature Tensor on a Manifold"
