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 hi_class . 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 8 \% 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 .