Abstract We study in $\mathbb {R}^{3+1}$ a system of nonlinearly coupled Klein-Gordon equations under the null condition, with (possibly vanishing) mass varying in the interval $[0, 1]$ . Our goal is three-fold, which extends the results in the earlier work of [5, 3]: 1) we want to establish the global well-posedness result to the system that is uniform in terms of the mass parameter (i.e., the smallness of the initial data is independent of the mass parameter); 2) we want to obtain a unified pointwise decay result for the solution to the system, in the sense that the solution decays more like a wave component (independent of the mass parameter) in a certain range of time, while the solution decays as a Klein-Gordon component with a factor depending on the mass parameter in the other part of the time range; 3) the solution to the Klein-Gordon system converges to the solution to the corresponding wave system in a certain sense when the mass parameter goes to 0. In order to achieve these goals, we will rely on both the flat and hyperboloidal foliation of the spacetime and prove a mass-independent $L^2$ –type energy estimate for the Klein-Gordon equations with possibly vanishing mass. In addition, the case of the Klein-Gordon equations with certain restricted large data is discussed.