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.