  We present a method to introduce relativistic corrections including linear
dark energy perturbations in Horndeski theory into Newtonian simulations based
on the N-body gauge approach. We assume that standard matter species (cold dark
matter, baryons, photons and neutrinos) are only gravitationally-coupled with
the scalar field and we then use the fact that one can include modified gravity
effects as an effective dark energy fluid in the total energy-momentum tensor.
In order to compute the scalar field perturbations, as well as the cosmological
background and metric perturbations, we use the Einstein-Boltzmann code
\hiclass. As an example, we study the impact of relativistic corrections on the
matter power spectrum in k-essence, a subclass of Horndeski theory, including
the effects of massless and massive neutrinos. For massive neutrinos with $\sum
m_{\nu} = 0.1$ eV, the corrections due to relativistic species (photons,
neutrinos and dark energy) can introduce a maximum deviation of approximately
$7\%$ to the power spectrum at $k \sim 10^{-3} \ \textrm{Mpc}^{-1}$ at $z=0$,
for a scalar field with sound speed $c_{s}^{2}\sim 0.013$ during matter
domination epoch. Our formalism makes it possible to test beyond $\Lambda$CDM
models probed by upcoming large-scale structure surveys on very large scales.
