The Lie derivative of the Christoffel symbol is
$$ \mathcal{L}_\xi \Gamma^k_{ij} = \nabla_i \nabla_j \xi^k - R^k_{ijl} \xi^l \,. $$
How can one prove that? And why does it make sense, because Christoffel symbols are functions? I know that the last question could be irrelevant, since the correct form of the LHS of the equation should be $(\mathcal{L}_\xi \Gamma)^k_{ij}$. But, I still cannot figure it out.
The formula gives the components of the Lie derivative of the connection as a whole, not the Lie derivative of each Christoffel symbol which is a function. Let's assume for a moment that the connection is a $(1,2)$ tensor and compute the Lie derivative formally. We have
$$ (\mathcal{L}_{\xi} \nabla)(X,Y) = \mathcal{L}_{\xi}(\nabla_X Y) - \nabla_{\mathcal{L}_{\xi} X} Y - \nabla_X (\mathcal{L}_{\xi} Y) = [\xi, \nabla_X Y] - \nabla_{[\xi, X]} Y - \nabla_X([\xi, Y]).$$
We have
$$ R^{\nabla}(X,Y) = \nabla_X \nabla_Y - \nabla_Y \nabla_X - \nabla_{[X,Y]} $$
and if $\nabla$ is torsion free also
$$ \nabla_{X} Y - \nabla_Y X = [X,Y]. $$
Using those formulas, we get
$$ (\mathcal{L}_{\xi}\nabla)(X,Y) = \\(\nabla_{\xi}(\nabla_X Y) - \nabla_{\nabla_X Y} \xi) + \\ (R^{\nabla}(\xi, X) Y - \nabla_{\xi}(\nabla_X Y) + \nabla_X (\nabla_{\xi} Y)) -\\ (\nabla_X(\nabla_{\xi} Y) - \nabla_X(\nabla_Y \xi)) = \\ R^{\nabla} (\xi, X) Y + \nabla_X (\nabla_Y \xi) - \nabla_{\nabla_{X} Y} \xi = R^{\nabla}(\xi, X) Y + (\nabla^2 \xi)(X,Y). $$
Writing the identity above in a local frame gets your equation with a plus sign between the two terms (this makes sense if your curvature tensor has the opposite sign convention).
To justify the derivation above, let's discuss how to define the Lie derivative of a connection. While a connection is not a tensor, the space of all connections form an affine space as the difference between two connections is a tensor. Given a diffeomorphism $\varphi \colon M \rightarrow M$ and a connection $\nabla$ on $TM$, we can get a new connection by the formula $$\nabla'_X Y = \varphi^{*}(\nabla_{\varphi_{*}(X)} (\varphi_{*}(Y))). $$ This is very similar to how a diffeomorphism acts on a $(1,2)$ tensor. Then we can mimick the definition of a the Lie derivative of a $(1,2)$ tensor and define
$$ (\mathcal{L}_{\xi}(\nabla)(p))(X,Y)|_{p} = \lim_{t \to 0} \frac{\left( \varphi_t^{*}(\nabla_{{\varphi_t}_{*}(X)}({\varphi_t}_{*}(Y))) - \nabla_X Y \right)|_{p}}{t}$$
where $\varphi_t$ is the local flow of $\xi$. Since the difference of two connections is a tensor, the expression above is well-defined and the result is a $(1,2)$-tensor which we can call the Lie derivative of $\nabla$.
BTW, a diffeomorphism $\varphi$ which preserves $\nabla$ in the sense that $\nabla_X Y = \varphi^{*}(\nabla_{\varphi_{*}(X)} (\varphi_{*}(Y)))$ is called an affine transformation and a vector field $\xi$ for which $\mathcal{L}_{\xi}(\nabla) = 0$ is sometimes called an affine vector field or an infinitesimal affine transformation because the flow of $\xi$ consists of affine transformations. This is analogous to Killing vector fields which are the infinitesimal generators of isometries and satisfy $\mathcal{L}_{\xi} g = 0$. For more information, you can read the section about affine transformations in the first volume of "Foundations of Differential Geometry" by Kobayashi and Nomizu.
