Recently, a rigorous renormalization theory for various scalar statistics has been developed for special modes of random advection diffusion involving random shear layer velocity fields with long-range spatiotemporal correlations. New random shearing direction models for isotropic turbulent diffusion are introduced here. In these models the velocity field has the spatial second-order statistics of an arbitrary prescribed stationary incompressible isotropic random field including long-range spatial correlations with infrared divergence, but the temporal correlations have finite range. The explicit theory of renormalization for the mean and second-order statistics is developed here. With epsilon the spectral parameter, for - infinity < epsilon < 4 and measuring the strength of the infrared divergence of the spatial spectrum, the scalar mean statistics rigorously exhibit a phase transition from mean-field behavior for epsilon < 2 to anomalous behavior for epsilon with 2 < epsilon < 4 as conjectured earlier by Avellaneda and the author. The universal inertial range renormalization for the second-order scalar statistics exhibits a phase transition from a covariance with a Gaussian functional form for epsilon with epsilon < 2 to an explicit family with a non-Gaussian covariance for epsilon with 2 < epsilon < 4. These non-Gaussian distributions have tails that are broader than Gaussian as epsilon varies with 2 < epsilon < 4 and behave for large values like exp(-C(epsilon) Absolute value of X 4-epsilon), With C(epsilon) an explicit constant. Also, here the attractive general principle is formulated and proved that every steady, stationary, zero-mean, isotropic, incompressible Gaussian random velocity field is well approximated by a suitable superposition of random shear layers.