We establish uniform error estimates of a finite difference method for the Zakharov system (ZS) with a dimensionless parameter $\epsilon$, which is inversely proportional to the acoustic speed. In the subsonic limit regime, the solution presents highly oscillatory initial layers due to the wave operator or the incompatibility of the initial data. Specifically, the solution propagates waves with $O(\epsilon)$-wavelength in time at O($\epsilon^2$) and O(1) amplitudes for well-prepared and ill-prepared initial data, respectively. In this talk, we rigorously analyze the error estimates for ZS via an asymptotic consistent formulation. For well-prepared initial data, we obtain error bounds at $O(h^2+t^{4/3})$ with time step t and mesh size h. While for ill-prepared initial data, the order is $O(h^2+t)$. The error estimates are confirmed by the numerical results.