The short answer
Eight order means it approximates $u$ with nine points and the local truncation error will be $O(\Delta{}x^8)$.
Higher the order, better the approximation.
The eight order mentionned is the accuracy order of the finite difference, nothing to do with the derivation order.
Have a look at this.
The long answer
Let us note $n$ the order of derivation and $m$ the order of accuracy of the method. You are solving the Navier-Stokes equation which contains first ($n=1$) and second ($n=2$) spatial derivatives.
DNS of the fully compressible Navier-Stokes equations for reacting flows are performed using explicit eight order central differences for spatial derivatives.
This means that you are using central differences of accuracy order $m=8$ to compute the $n=1$ and $n=2$ order spatial derivatives. The greater the order $m$, the more precise the method because the order quantify the local truncation error. To compute an finite difference expression of accuracy order $m$ you need $m+1$ points. For example, let us consider that we want a forward finite difference of order 1 for the first spatial derivative : we need two points, $x_0$ and $x_0+\Delta{}x$. The Taylor expansion of the function $u$ at $x_0+\Delta{}x$ reads
$$u(x_0+\Delta{}x)=u(x_0)+\Delta{}xu_x(x_0)+O(\Delta{}x^2).$$
With a little bit of algebra we get
$$u_x(x_0)=\frac{u_0(x+\Delta{}x)-u(x_0)}{\Delta{}x}+\underbrace{O(\Delta{}x)}_{\begin{matrix}\mbox{truncation}\\\mbox{error}\end{matrix}}$$
and the resulting truncation error gives us the accuracy order of the finite difference. Here the local truncation error is proportional to the step size. To have a good accuracy it will therefore be necessary to have a very fine mesh size which is very expensive in calculation. While with a finite difference of order $m=2$ and with a truncation error $O(\Delta{}x^2)$, the mesh size may be larger for the same numerical error.
It exists some tools as this one to compute the finite difference coefficient of order $m$. If you click on "How does it work?" you will have an explanation of how to do it by yourself. You may also be interested in this Wikipedia article providing the coefficient for some high order finite differences. According to the latter, the $m=8$ central finite difference for $n=1$ spatial derivative of $u$ is given by
\begin{align*}
u_{x,i}&=\frac{3u_{i-4}-32u_{i-3}+168u_{i-2}-672u_{i-1}+672u_{i+1}-168u_{i+2}+32u_{i+3}-3u_{i+4}}{840\Delta{x}}\\&+O{\Delta{x}^8}
\end{align*}
with $u_i = u(x+i\Delta{x})$.
Edit: precision about explicit filtering by Johntra Volta
Johntra Volta added in the comment some interesting explanation about explicit filtering so I added it to my answer for more visibility.
Short answer:
Perform a convolution with a low pass filter.
Long answer:
Even if we use all the points in one direction to form FD stencil, (spectral differentiation) we may still not be able to prevent blow-up of simulation after longer time-integration. The reason is aliasing of unresolved spatial scales to the ones represented by the numerical grid. The build-up of energy due to aliasing starts from the highest end of resolved scales (appearance of "2h waves"). One remedy is low-pass filtering of the solution. You may use same web app to create kernel coefficients (Hint: put $n$=2).