We propose and rigourously analyze a multiscale time integrator Fourier pseudospectral (MTI-FP) method for the Dirac equation with a dimensionless parameter \(\varepsilon\in(0,1]\) which is inversely proportional to the speed of light. In the nonrelativistic limit regime, i.e. \(0<\varepsilon\ll 1\), the solution exhibits highly oscillatory propagating waves with wavelength \(O(\varepsilon^2)\) and \(O(1)\) in time and space, respectively. Due to the rapid temporal oscillation, it is quite challenging in designing and analyzing numerical methods with uniform error bounds in \(\varepsilon\in(0,1]\). We present the MTI-FP method based on properly adopting a multiscale decomposition of the solution of the Dirac equation and applying the exponential wave integrator with appropriate numerical quadratures. By a careful study of the error propagation and using the energy method, we establish two independent error estimates via two different mathematical approaches as \(h^{m_0}+\frac{\tau^2}{\varepsilon^2}\) and \(h^{m_0}+\tau^2+\varepsilon^2\), where \(h\) is the mesh size, \(\tau\) is the time step and \(m_0\) depends on the regularity of the solution. These two error bounds immediately imply that the MTI-FP method converges uniformly and optimally in space with exponential convergence rate if the solution is smooth, and uniformly in time with linear convergence rate at \(O(\tau)\) for all \(\varepsilon\in(0,1]\) and optimally with quadratic convergence rate at \(O(\tau^2)\) in the regimes when either \(\varepsilon=O(1)\) or \(0<\varepsilon\lesssim \tau\). Numerical results are reported to demonstrate that our error estimates are optimal and sharp. Finally, the MTI-FP method is applied to study numerically the convergence rates of the solution of the Dirac equation to those of its limiting models when \(\varepsilon\to0^+\).