A multiscale time integrator Fourier pseudospectral (MTI-FP) method is presented for discretizing the massive Klein-Gordon-Dirac (KGD) system which involves a small dimensionless parameter 0<ε≤1. In the nonrelativistic limit regime, the KGD system admits rapid oscillations in time as ε→0+. In addition, the nonlinear Yukawa interaction and the indefinite Dirac operator bring other significant difficulties. The main idea of the MTI-FP method is to construct a precise multiscale decomposition by the frequency (MDF) to the solution of the KGD system at each time step and then employ the Fourier pseudospectral discretization for the spatial derivatives followed with the exponential wave integrator (EWI) for the time marching. This approach is explicit, easy to implement and preforms significantly better than the classical methods in the literature. More specifically, we rigorously establish the uniform error bounds at O(τ+hm0−1) for all ε∈(0,1] and optimal quadratic temporal error bounds at O(τ2) in the ε=O(1) regime, where τ is the time step size, h is mesh size and m0 depends on the regularity of the solution. Extensive numerical results demonstrate that our error bounds are optimal and sharp. Finally, we apply the MTI-FP method to numerically study the nonrelativistic limit behaviors of the KGD system when ε→0+.