Consider a metric manifold $(M,g_{\mu\nu})$. Consider also a connection $\Gamma^{\mu}_{\alpha\beta}$ defined on it, in turn defining a covariant derivative $$ \nabla_\mu e_{\nu} = \Gamma_{\mu\nu}^\rho e_\rho, $$ without making any assumption on the symmetry of its lower indices (i.e. without the hypothesis of vanishing torsion).
Can we infer from the variation of the Einstein-Hilbert action $$ S_{EH}[g_{\mu\nu}, \Gamma^{\rho}_{\alpha\beta}] = \int d^D\!x \sqrt{-g} R_{\mu\nu} g^{\mu\nu}, $$ with respect to $\Gamma^{\rho}_{\alpha\beta}$, that the connection should satisfy both $$ \Gamma^{\rho}_{\mu\nu}=\Gamma^{\rho}_{\nu\mu}\qquad \text{and}\qquad \nabla_{\alpha}g_{\mu\nu}=0? $$ I tried using the Palatini formula $$ \delta R_{\mu\nu} = \nabla_\nu \delta \Gamma^{\sigma}_{\sigma \mu}- \nabla_\sigma \delta\Gamma^{\sigma}_{\nu\mu} + 2 S^{\rho}_{\sigma\nu}\delta\Gamma^{\sigma}_{\rho \mu}, $$ where $S^\alpha_{\mu\nu}=\Gamma^\alpha_{[\mu\nu]}$ is torsion, which gives $$ \delta S_{EH} = \int d^D\!x \sqrt{-g} \left[ \nabla_\nu \delta \Gamma^{\sigma}_{\sigma \mu}- \nabla_\sigma \delta\Gamma^{\sigma}_{\nu\mu} + 2 S^{\rho}_{\sigma\nu}\delta\Gamma^{\sigma}_{\rho \mu} \right] g^{\mu\nu} $$ and integrating by parts $$ \delta S_{EH} = \int d^D\!x \delta \Gamma^{\sigma}_{\nu\mu} \left[ \nabla_{\sigma} \left( \sqrt{-g} g^{\mu\nu} \right)- \delta^\nu_\sigma \nabla_\rho\left(\sqrt{-g}g^{\sigma\mu}\right) + 2S^{\nu}_{\alpha\sigma} g^{\alpha\mu} \right]. $$ Now I can use $$\nabla\sqrt{-g} = -\frac{1}{2}\sqrt{-g} g_{\alpha\beta}\nabla_\sigma g^{\alpha\beta},$$ but I seem to get nowhere.
