Post-Newtonian celestial dynamics is a relativistic theory of motion of massive bodies and test particles under the influence of relatively weak gravitational forces. Standard approach for development of this theory relies upon the key concept of the isolated astronomical system supplemented by the assumption that the background space-time is flat. The standard post-Newtonian theory of motion was instrumental in explanation of the existing experimental data on binary pulsars, satellite and lunar laser ranging, and in building precise ephemerides of planets in the solar system. Recent studies of the formation of large-scale structure in our universe indicate that the standard post-Newtonian mechanics fails to describe more subtle dynamical effects in motion of the bodies comprising the astronomical systems of larger size -galaxies and clusters of galaxies -where the Riemann curvature of the expanding FLRW universe interacts with the local gravitational field of the astronomical system and, as such, can not be ignored.The present paper outlines theoretical principles of the post-Newtonian mechanics in the expanding universe. It is based upon the gauge-invariant theory of the Lagrangian perturbations of cosmological manifold caused by an isolated astronomical N-body system (the solar system, a binary star, a galaxy, a cluster of galaxies). We postulate that the geometric properties of the background manifold are described by a homogeneous and isotropic Friedman-Lemaître-Robertson-Walker (FLRW) metric governed by two primary components -the dark matter and the dark energy. The dark matter is treated as an ideal fluid with the Lagrangian taken in the form of pressure along with the scalar Clebsch potential as a dynamic variable. The dark energy is associated with a single scalar field with a potential which is hold unspecified as long as the theory permits. Both the Lagrangians of the dark matter and the scalar field are formulated in terms of the field variables which play a role of generalized coordinates in the Lagrangian formalism. It allows us to implement the powerful methods of variational calculus to derive the gauge-invariant field equations of the post-Newtonian celestial mechanics of an isolated astronomical system in an expanding universe. These equations generalize the field equations of the post-Newtonian theory in asymptotically-flat spacetime by taking into account the cosmological effects explicitly and in a self-consistent manner without assuming the principle of liner superposition of the fields or a vacuole model of the isolated system, etc. The field equations for matter dynamic variables and gravitational field perturbations are coupled in the most general case of arbitrary equation of state of matter of the background universe.We introduce a new cosmological gauge which generalizes the de Donder (harmonic) gauge of the post-Newtonian theory in asymptotically flat spacetime. This gauge significantly simplifies the gravitational field equations and allows to find out the approximations where the ...