I am deriving the Einstein equation using the Einstein-Hilbert action: 
It is obvious that the variation in the Riemann Tensor is calculated from a variational product rule. What is not obvious to me is why variations obey this rule, and I'll like an explanation.
The easiest way to go about this is to think of variation in terms of the background field method. Here, you write the metric as some background plus perturbation:
\begin{equation} g_{\mu \nu} = \bar{g}_{\mu \nu} + h_{\mu \nu} \end{equation}
and then you expand every quantity that depends on $g_{\mu \nu}$ as a power series in $h$:
$$ \Gamma = \bar{\Gamma} + \Gamma^{(1)} + \Gamma^{(2)} \cdots \\ R = \bar{R} + R^{(1)} + R^{(2)} \cdots $$
where the superscript $^{(n)}$ labels the $\mathcal{O} (h^n)$ perturbation. You can then easily see that
$$ \Gamma \Gamma = (\bar{\Gamma} + \Gamma^{(1)} + \Gamma^{(2)} \cdots)(\bar{\Gamma} + \Gamma^{(1)} + \Gamma^{(2)} \cdots) $$
Extracting the $\mathcal{O}(h)$ variation, like in your problem, we automatically get both the contributions:
$$ (\Gamma \Gamma)^{(1)} = \bar{\Gamma} \Gamma^{(1)} + \Gamma^{(1)} \bar{\Gamma} $$
It's a fundamental property of infinitesimal variations $\delta$ that they are linear derivations, i.e. they obey Leibniz product rule. Short of a rigorous definition of infinitesimals, the proof is basically to consider variations as 1-parameter families of functions, and then differentiate wrt. the parameter, cf. above comment by Bence Racsko.
