A two-level fully discrete finite element variational multiscale method based on two local Gauss integrations for the unsteady Navier-Stokes equations is presented and studied, where conforming finite element pairs are used for the spatial discretization and a three-point difference formula is employed for the temporal discretization. At each time step of this method, a stabilized nonlinear Navier-Stokes system is first solved on a coarse grid, and then a stabilized linear problem is solved on a fine grid to correct the coarse grid solution, where the numerical forms of the Navier-Stokes equations both on coarse and fine meshes are stabilized by a stabilization term defined by the difference between a consistent and an under-integrated matrix of the gradient of velocity interpolant. Error bound of the approximate velocity is derived. Numerical results are also given to verify the theoretical predictions and demonstrate the effectiveness of the proposed numerical method.