  The N-body gauge allows the introduction of relativistic effects in Newtonian
cosmological simulations. Here we extend this framework to general Horndeski
gravity theories, and investigate the relativistic effects that the scalar
field introduces in the matter power spectrum at intermediate and large scales.
In particular, we show that the kineticity function at these scales enhances
the amplitude of the signal of contributions coming from the extra degree of
freedom. Using the Quasi-Static Approximation (QSA), we separate modified
gravity effects into two parts: one that only affects small-scale physics, and
one that is due to relativistic effects. This allows our formalism to be
readily implemented in modified gravity N-body codes in a straightforward
manner, e.g., relativistic effects can be included as an additional linear
density field in simulations. We identify the emergence of gravity acoustic
oscillations (GAOs) in the matter power spectrum at large scales, $k \sim
10^{-3}-10^{-2}$ Mpc$^{-1}$. GAO features have a purely relativistic origin,
coming from the dynamical nature of the scalar field. GAOs may be enhanced to
detectable levels by the rapid evolution of the dark energy sound horizon in
certain modified gravity models and can be seen as a new test of gravity at
scales probed by future galaxy and intensity-mapping surveys.
