A numerical study of fundamental aspects of turbulent mixing has been performed,with emphasis on the behavior of passive scalars of low molecular diffusivity (high Schmidt number Sc). Direct Numerical Simulation is used to simulate incompressible, stationary and isotropic turbulence carried out at high grid resolution. Data analyses are carried out by separate parallel codes using up to 1024^3 grid points for Taylor-scale Reynolds number (R_lambda) up to 390 and Sc up to 1024.Schmidt number of order 1000 is simulated using a double-precision parallel code in a turbulent flow at a low Reynolds number of R_lambda 8 to reduce computational cost to achievable level. The results on the scalar spectrum at high Schmidt numbers appear to have a k^{-1} scaling range. In the presence of a uniform mean scalar gradient, statistics of scalar gradients are observed to deviate substantially from Kolmogorov's hypothesis of local isotropy, with a skewness factor remaining at order unity as the Reynolds number increases. However, this skewness decreases with Schmidt number suggesting that local isotropy for scalars at high Schmidt number is a better approximation. Intermittency exponents manifested by three types of two-point statistics of energy and scalar dissipation, i.e., the two-point correlator (chi(x)chi(x+r)), the second-order moment of local scalar dissipation (chi_r^2) and the variance of the logarithmic local scalar dissipation sigma^2_{lnchi_r} are discussed. Several basic issues in differential diffusion between two scalars of different molecular diffusivities transported by the same turbule nt flow, the physical process of scalar spectral transfer and subgrid-scale transfer are also briefly addressed.