Post-Newtonian celestial dynamics is a relativistic theory of motion of massive bodies and test particles under the influence of relatively weak gravitational forces. The standard approach for development of this theory relies upon the key concept of the isolated astronomical system supplemented by the assumption that the background spacetime is flat. The standard post-Newtonian theory of motion was instrumental in the 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 structures 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 Friedmann-Lema\^{\i}tre-Robertson-Walker universe interacts with the local gravitational field of the astronomical system and, as such, cannot 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, and a cluster of galaxies). We postulate that the geometric properties of the background manifold are described by a homogeneous and isotropic Friedmann-Lema\^{\i}tre-Robertson-Walker 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 an 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 one to find out the approximations where the field equations can be fully decoupled and solved analytically. The residual gauge freedom is explored and the residual gauge transformations are formulated in the form of the wave equations for the gauge functions. We demonstrate how the cosmological effects interfere with the local system and affect the local distribution of matter of the isolated system and its orbital dynamics. Finally, we worked out the precise mathematical definition of the Newtonian limit for an isolated system residing on the cosmological manifold. The results of the present paper can be useful in the Solar System for calculating more precise ephemerides of the Solar System bodies on extremely long time intervals, in galactic astronomy to study the dynamics of clusters of galaxies, and in gravitational wave astronomy for discussing the impact of cosmology on generation and propagation of gravitational waves emitted by coalescing binaries and/or merging galactic nuclei.