I came to this post many years later, since I too was concerned about the absence of Riemann's formula in most texts, lengthy treatment in others, or reliance on more advanced techniques like Jacobi fields. I include here a direct concise proof which I think would be well suited for beginning students. We want to show that
$$
g_{ij}(x)=\delta_{ij}+\frac{1}{3} R_{k ij\ell}(o)x_kx_\ell+\mathcal{O}(|x|^3),
$$
where $x=(x_1,\dots, x_n)$ are normal coordinates centered at a point $o$ in a Riemannian manifold $M$. By Taylor's theorem, we need to check that
\begin{eqnarray}
\tag{1} g_{ij}(o)&=&\delta_{ij}\label{eq:1},\\
\tag{2} g_{ij,k}(o)&=&0\label{eq:2},\\
\tag{3} g_{ij,k\ell}(o)&=&\frac{1}{3}\big(R_{kij\ell }(o)+R_{\ell ij k}(o)\big)\label{eq:3},
\end{eqnarray}
where we use the notation $f_{,i}:=\partial_i f$, and $f_{,ij}:=\partial^2_{ij} f$.
\eqref{eq:1} follows immediately from the construction of normal coordinates, and \eqref{eq:2} is not difficult to establish either. \eqref{eq:3}, which is the heart of the matter, requires a bit more work.
In most sources, like Spivak, vol II or the relatively recent book by Jost, which gives a very comprehensive treatment of Riemann's lecture, \eqref{eq:3} is established via symmetry properties of $g_{ij,k\ell}$ which involve long computations; although 2011 Lecture Notes of John Douglas Moore gives a very nice and efficient proof of them. Instead I will demonstrate \eqref{eq:3} via a Bianchi type cyclic identity for derivatives of Christoffel symbols, which follows quickly from the geodesic equation (as the OP had mentioned in the "Update/Scholium" above). The only place I have seen this approach is in 2013 Lecture Notes of Christian Bär.
Proofs of (1) and (2)
All indices here range from $1$ to $n$, and any term which involves repeated indices stands for a sum over that index.
Let $\exp_o\colon T_o M\to M$ be the exponential map, and $U\subset M$ be a ball centered at $o$ such that $\exp_o\colon \exp_o^{-1}(U)\to U$ is a diffeomorphism. Let $e_i$ be an orthonormal basis for $T_o M$, i.e.,
$
g(e_i, e_j)=\delta_{ij}.
$
Then the normal coordinates $x_i\colon U\to \mathbf{R}$ (with respect to $e_i$) are given by
$$
\exp_o(x_i(p) e_i)=p.
$$
The mapping $x\colon U\to\mathbf{R}^n$, given by $x:=(x_1,\dots, x_n)$ identifies $U$ with a ball centered at the origin in $\mathbf{R}^n$, which we again denote by $U$. Let $E_i(x):=\partial_i|_x$ be the coordinate vector fields on $U$. Then $g_{ij}\colon U\to\mathbf{R}$ are given by
$$
g_{ij}(x):=g(E_i(x),E_j(x)).
$$
Since $E_i(o)=e_i$ we immediately obtain \eqref{eq:1}. To see \eqref{eq:2} note that, by the definition of normal coordinates $x_i$, the geodesics in $U$ passing through $o$ and another point $x$ of $U$ are given by $\gamma(t):=tx$. Since $\gamma$ is a geodesic,
$$
\gamma_k''(t)+\Gamma_{ij}^k(\gamma(t))\gamma_i'(t)\gamma_j'(t)=0,
$$
which yields
\begin{equation}\label{eq:4}
\tag{4}\Gamma_{ij}^k(t x)x_i x_j=0.
\end{equation}
Setting $t=0$, observing that $x_i$, $x_j$ may assume any values, and recalling that $\Gamma_{ij}^k=\Gamma_{ji}^k$, we obtain
\begin{equation*}\label{eq:Gamma0}
\Gamma_{ij}^k(o)=0.
\end{equation*}
Since
$
\Gamma_{ij}^k=\frac{1}{2} g^{k\ell}(g_{\ell i,j}+g_{\ell j,i}-g_{ij,\ell}),
$
and $g^{ij}(o)=\delta_{ij}$, we have
\begin{equation*}\label{eq:8}
0=\Gamma_{ij}^k(o)=\frac{1}{2} \big(g_{k i,j}(o)+g_{k j,i}(o)-g_{ij,k}(o)\big).
\end{equation*}
Adding the above equation to itself, after a cyclic permutation of indices yields \eqref{eq:2}.
Proof of (3)
Since
$
\nabla_{E_i}E_j=\Gamma_{ij}^kE_k,
$
$$
g_{ij,k}=g(\nabla_{E_k}E_i, E_j)+g(E_i, \nabla_{E_k}E_j)=\Gamma_{k i}^\ell g_{\ell j}+\Gamma_{k j}^\ell g_{i\ell}.
$$
Differentiating again, and using \eqref{eq:2}, yields
\begin{equation}\label{eq:5}
\tag{5}g_{ij,k\ell}(o)=\Gamma_{k i,\ell}^j(o)+\Gamma_{k j,\ell}^i(o).
\end{equation}
Next we differentiate \eqref{eq:4} at $t=0$ to obtain
$$
\Gamma_{ij,\ell}^k(o)x_i x_jx_\ell=0,
$$
a homogeneous polynomial of degree $3$ which vanishes identically. The coefficient of each term $x_i x_jx_\ell$ is the sum of all $6$ permutations of lower indices of $\Gamma_{ij,\ell}^k(o)$. Since $\Gamma_{ij,\ell}^k=\Gamma_{ji,\ell}^k$, we obtain
\begin{equation*}\label{eq:cyclic}
\Gamma_{ij,\ell}^k(o)+\Gamma_{j\ell,i}^k(o)+\Gamma_{\ell i,j}^k(o)=0.
\end{equation*}
Now note that, since
$
R_{ijk}^\ell=\Gamma_{ik,j}^\ell-\Gamma_{jk,i}^\ell+\Gamma_{ik}^p\Gamma_{pj}^\ell-\Gamma_{jk}^p\Gamma_{ip}^\ell,
$
\eqref{eq:2} yields that
$$
R_{ijk\ell}(o)=\Gamma_{ik,j}^\ell(o)-\Gamma_{jk,i}^\ell(o).
$$
Here we have also used the fact that $R_{ijk}^\ell(o)=R_{ijk\ell}(o)$ due to \eqref{eq:1}. The last two displayed equations yield
\begin{eqnarray*}
R_{ik\ell j}(o)+R_{i\ell kj}(o)
= \Gamma_{i\ell,k}^j(o)+\Gamma_{ik,\ell}^j(o)-2\Gamma_{\ell k,i}^j(o)
=-3\Gamma_{\ell k,i}^j(o).
\end{eqnarray*}
The last equality together with \eqref{eq:5} and symmetries of $R$ now yields
\begin{eqnarray*}
g_{ij,k\ell}(o)
&=&-\frac{1}{3}\big(R_{\ell ik j}(o)+\require{cancel}\cancel{R_{\ell kij}(o)}+R_{\ell jk i}(o)+\require{cancel}\cancel{R_{\ell kji}(o)\big)}\\
&=&\frac{1}{3}\big(R_{\ell i jk}(o)+R_{kij\ell}(o)\big),
\end{eqnarray*}
as desired.
