Different efficient and accurate numerical methods have recently been proposed and analyzed for the nonlinear Klein-Gordon equation (NKGE) with a dimensionless parameter ε ∈ (0, 1], which is inversely proportional to the speed of light. In the nonrelativestic limit regime, i.e. 0 < ε 1, the solution of the NKGE propagates waves with wavelength at O(1) and O(ε 2 ) in space and time, respectively, which brings significantly numerical burdens in designing numerical methods. We compare systematically spatial/temporal efficiency and accuracy as well as ε-resolution (or ε-scalability) of different numerical methods including finite difference time domain methods, time-splitting method, exponential wave integrator, limit integrator, multiscale time integrator, two-scale formulation method and iterative exponential integrator. Finally, we adopt the multiscale time integrator to study the convergence rates from the NKGE to its limiting models when ε → 0 + .