Does one need to take into account finite gravity speed in N-body simulations? Most small-scale N-body simulations (e.g., planetary systems, solar system, stellar clusters, ...) use  classical Newtonian gravity. 

Most large-scale N-body simulations (e.g., galaxy clusters, ...) use general relativity.

My question: if I were to make a large-scale N-body simulation using Newtonian gravity, would it be enough to "correct" that to first-order by modifying the equations of motion to take into account the fact that gravitational fields propagate at the speed of light?
If you're asking whether it's sufficient to use a retarded (time-delayed) positions to calculate gravitational forces, then no, that would be much worse than Newtonian gravity. For example, that would predict that the Earth should spiral into the Sun on the order of about 400 years. See also answers to 


Besides retarded gravitation, anything else to worry about when calculating MU69's orbit from scratch?
Is the influence of gravity instantaneous?





  Most small-scale N-body simulations (e.g., planetary systems, solar system, stellar clusters, ...) use classical Newtonian gravity.


This is true, but some simulations use relativistic corrections.

The post-Newtonian expansion works in order $\epsilon\sim v^2\sim U$, where $U$ is (the negative of) Newtonian potential, and including the first-order post-Newtonian correction to the N-body equations of motion forms the Einstein–Infeld–Hoffmann equations. Because they contain $3$-body terms, the simulation would now be $\mathcal{O}(N^3)$ instead of the Newtonian simulation $\mathcal{O}(N^2)$.

The post-Newtonian three-body equations are known to 2PN, while the two-body problem to 3.5PN. Effects beyond the first order can be important in e.g. binary simulations because gravitational radiation losses are seen at 2.5PN order. However, they are not translatable to an N-body problem.

A technique of making the EIH equations more tractable by making a second expansion that depends on the problem at hand can be found in Will (2014), which is applicable when most the relativistic dynamics are dominated by a few contribution (e.g., a supermassive central black hole in galaxy simulations).


  Most large-scale N-body simulations (e.g., galaxy clusters, ...) use general relativity.


As Walter says, most of them use Newtonian gravity on a expanding background. Not necessarily directly computer particle-particle Newtonian force, but perhaps modeling Poisson's equation through mesh methods, or perhaps some combination thereof (e.g., direct force calculations for close particles, mesh methods for larger scales).


  My question: if I were to make a large-scale N-body simulation using Newtonian gravity, would it be enough to "correct" that to first-order by modifying the equations of motion to take into account the fact that gravitational fields propagate at the speed of light?


If you discount cosmological expansion, it will simply be wrong. However, first-order corrections to the Newtonian gravity part of simulations do exist in the post-Friedmann formalism. I'm not really familiar with PFF, but according to Bruni et al. (2014), the leading-order correction to Newtonian gravity introduces a vector potential in the $g_{0i}$ term in the metric that can produce weak lensing effects but doesn't affect matter dynamics. Which is entirely sensible because that kind of gravitomagnetic effects should be suppressed by a further $v/c$ term for matter, just like magnetism.

For references on how low order GR corrections to Newtonian gravity simulations are often made, see answers to How to calculate the planets and moons beyond Newtons's gravitational force?



References:


Will, C. M. "Incorporating post-Newtonian effects in N-body dynamics". Phys. Rev. D 89, 044043 (2014) [arXiv:1312.1289]
Bruni, M., Thomas, D. B., Wands, D. "Computing general-relativistic effects from Newtonian N-body simulations: Frame dragging in the post-Friedmann approach." Phys. Rev. D 89, 044010 (2014) [arXiv:1306.1562]