Not yet reviewed by Pith; the record is open.

We provide two derivations of the baryonic equations that can be straightforwardly implemented in existing Einstein--Boltzmann solvers. One of the derivations begins with an action principle, while the other exploits the conservation of the stress-energy tensor. While our result is manifestly covariant and satisfies the Bianchi identities, we point out that this is not the case for the implementation of the seminal work by Ma and Bertschinger and in the existing Boltzmann codes. We also study the tight coupling approximation up to the second order without choosing any gauge using the covariant full baryon equations. We implement the improved baryon equations in a Boltzmann code and investigate the change in the estimate of cosmological parameters by performing an MCMC analysis. With the covariantly correct baryon equations of motion, we find 1% deviation for the best fit values of the cosmological parameters that should be taken into account. While in this paper, we study the Lambda-CDM model only, our baryon equations can be easily implemented in other models and various modified gravity theories.