Aiming at discriminating different gravitational potential models of the Milky Way, we perform tests based on the kinematic data powered by the Gaia DR2 astrometry over a large range of (R, z) locations. Invoking the complete form of Jeans equations that admit three integrals of motion, we use the independent R- and z-directional equations as two discriminators (T-R and T-z). We apply the formula for spatial distributions of radial and vertical velocity dispersions proposed by Binney et al., and successfully extend it to azimuthal components, sigma(theta)(R, z) and V-theta(R, z); the analytic form avoids the numerical artifacts caused by numerical differentiation in Jeans-equations calculation given the limited spatial resolutions of observations, and more importantly reduces the impact of kinematic substructures in the Galactic disc. It turns out that whereas the current kinematic data are able to reject Moffat's Modified Gravity (let alone the Newtonian baryon-only model), Milgrom's MOND is still not rejected. In fact, both the carefully calibrated fiducial model invoking a spherical dark matter (DM) halo and MOND are equally consistent with the data at almost all spatial locations (except that probably both have respective problems at low-|z| locations), no matter which tracer population or which meaningful density profile is used. Since there is no free parameter at all in the quasi-linear MOND model we use, and the baryonic parameters are actually fine-tuned in the DM context, such an effective equivalence is surprising, and might be calling forth a transcending synthesis of the two paradigms.