Varying the Ricci tensor with respect to the metric $g^{\mu\nu}$, one would get
$ \delta R_{\mu\nu} = \delta(g^{\rho\sigma} R_{\rho\mu\sigma\nu}) = g^{\rho\sigma}\delta R_{\rho\mu\sigma\nu} + \delta g^{\rho\sigma} R_{\rho\mu\sigma\nu}, $
but in all my references I found that $\delta R_{\mu\nu} = g^{\rho\sigma}\delta R_{\rho\mu\sigma\nu}$ (including Wikipedia), which implies that the last term of the above equation is identically zero. Why is that so?
Consider a metric $\widetilde{g}_{\mu\nu}=g_{\mu\nu}+\delta g_{\mu\nu}$ and its inverse metric $\widetilde{g}^{\mu\nu}=g^{\mu\nu}-\delta g^{\mu\nu}$. (Fittingly I've just explained the negative sign here yesterday.) We have: \begin{align*} R_{\mu\nu}[\widetilde{g}] &=\widetilde{g}^{\kappa\lambda} R_{\kappa\mu\lambda\nu}[\widetilde{g}] =(g^{\kappa\lambda}-\delta g^{\kappa\lambda}) (R_{\kappa\mu\lambda\nu}[g]+\delta R_{\kappa\mu\lambda\nu}) \\ &=R_{\mu\nu}[g] +\underbrace{g^{\kappa\lambda}\delta R_{\kappa\mu\lambda\nu} -\delta g^{\kappa\lambda}R_{\kappa\mu\lambda\nu}[g]}_{=\delta R_{\mu\nu}}+\mathcal{O}(\delta g^2). \end{align*} I therefore don't think this formula is correct in a stricter sense, somebody probably approximated $\widetilde{g}^{\kappa\lambda}\approx g^{\kappa\lambda}$ directly in the first step and therefore left out a term in first order. Because of this I'd also use a different formula for the pertubation of the curvature tensor. Taking the formula for the Christoffel symbols and putting the perturbed metric in, you get: \begin{equation} \Gamma_{\mu\nu}^\kappa[\widetilde{g}] =\Gamma_{\mu\nu}^\kappa[g] +\underbrace{\frac{1}{2}g^{\kappa\lambda}\left( \nabla_\mu\delta g_{\lambda\nu} +\nabla_\nu\delta g_{\lambda\mu} -\nabla_\lambda\delta g_{\mu\nu}\right)}_{=\delta\Gamma_{\mu\nu}^\kappa} +\mathcal{O}(\delta g^2). \end{equation} This calculation is not short, but pretty straight forward. You immediately see, that $\delta\Gamma_{\mu\nu}^\kappa$ is a tensor, as it is a sum of covariant differentiations of a tensor. Taking the formula for the Riemann curvature tensor and putting this equation in, you get: \begin{equation} R_{\lambda\mu\nu}^\rho[\widetilde{g}] =R_{\lambda\mu\nu}^\rho[g] +\underbrace{\partial_\mu\delta\Gamma_{\lambda\nu}^\rho -\partial_\nu\delta\Gamma_{\lambda\mu}^\rho +\delta\Gamma_{\mu\sigma}^\rho\Gamma_{\lambda\nu}^\sigma[g] +\Gamma_{\mu\sigma}^\rho[g]\delta\Gamma_{\lambda\nu}^\sigma -\delta\Gamma_{\nu\sigma}^\rho\Gamma_{\lambda\mu}^\sigma[g] -\Gamma_{\nu\sigma}^\rho[g]\delta\Gamma_{\lambda\mu}^\sigma}_{=\delta R_{\lambda\mu\nu}^\rho} +\mathcal{O}(\delta g^2). \end{equation} Notice, that you can shorten this using the covariant differentiation: \begin{equation} \delta R_{\lambda\mu\nu}^\rho =\nabla_\mu\delta\Gamma_{\lambda\nu}^\rho -\nabla_\nu\delta\Gamma_{\lambda\mu}^\rho. \end{equation} You again immediately see, that $\delta R_{\lambda\mu\nu}^\rho$ is a tensor, as it is the sum of covariant differentiations of a tensor. Contracting to $\delta R_{\lambda\nu}$ gives you the simple Palatini identity.
This is for whoever is interested in the solution. The second term in the equation of my original post is not zero as it is. But, from $g_{\mu\alpha}g^{\alpha\nu}=\delta_\mu^\nu$, one can prove that $\delta g_\rho^\sigma = - g_\alpha^\sigma g_\rho^\beta \delta g_\beta^\alpha$. Then the variation of the Ricci tensor gives
$ \delta R_{\mu\nu} = \delta(g_\rho^\sigma R^\rho_{\;\mu\sigma\nu})= g^\sigma_\rho\delta R^\rho_{\;\mu\sigma\nu} + \delta g_\rho^\sigma R^\rho_{\;\mu\sigma\nu}, $
where
$ \delta g_\rho^\sigma R^\rho_{\;\mu\sigma\nu} = - g_\alpha^\sigma g_\rho^\beta \delta g_\beta^\alpha R^\rho_{\;\mu\sigma\nu} = -\delta g_\beta^\alpha R^\beta_{\;\mu\alpha\nu}, $
which vanishes as $\rho,\sigma,\alpha,\beta$ are all dummy indices. Hence,
$ \delta R_{\mu\nu} = g^\sigma_\rho\delta R^\rho_{\;\mu\sigma\nu}. $
The variation of the Ricci tensor $$ \delta R_{\mu\nu}~=~\delta(g^{\rho\sigma} R_{\rho\mu\sigma\nu}) = g^{\rho\sigma}\delta R_{\rho\mu\sigma\nu} + \delta g^{\rho\sigma} R_{\rho\mu\sigma\nu}~=~g^{\rho\sigma}\delta R_{\rho\mu\sigma\nu}~+~\delta s\left(\frac{D}{ds} g^{\rho\sigma}\right)R_{\rho\mu\sigma\nu}. $$ The last term is zero by the covariant constancy of the metric. This is seen with $$ \frac{D}{ds} g^{\rho\sigma}~=~\frac{d}{ds} g^{\rho\sigma}~+~\left(\Gamma^\rho_{\alpha\beta}g^{\alpha\sigma}~+~\Gamma^\sigma_{\alpha\beta}g^{\alpha\rho}\right)\frac{dx^\beta}{ds}. $$ This is why the $\delta g^{\rho\sigma} R_{\rho\mu\sigma\nu}~=~0$
