  The probability density function (PDF) of the logarithmic density contrast,
$s=\ln (\rho/\rho_0)$, with gas density $\rho$ and mean density $\rho_0$, for
hydrodynamical supersonic turbulence is well-known to have significant
non-Gaussian (intermittent) features that monotonically increase with the
turbulent Mach number, $\mathcal{M}$. By studying the mass- and volume-weighted
$s$-PDF for an ensemble of 36 sub-to-trans-Alf\'venic mean-field, supersonic,
isothermal turbulence simulations with different modes of driving, relevant to
molecular gas in the cool interstellar medium, we show that a more intricate
picture emerges for the non-Gaussian nature of $s$. Using four independent
measures of the non-Gaussian components, we find hydrodynamical-like structure
in the highly magnetised plasma for $\mathcal{M} \lesssim 4$. However, for
$\mathcal{M} \gtrsim 4$, the non-Gaussian signatures disappear, leaving
approximately Gaussian $s$-statistics -- exactly the opposite of hydrodynamical
turbulence in the high-$\mathcal{M}$ limit. We also find that the non-Gaussian
components of the PDF increase monotonically with more compressive driving
modes. To understand the $\mathcal{M} \lesssim 4$ non-Gaussian features we use
one-dimensional (1D) pencil beams to explore the dynamics along and across the
large-scale magnetic field, $\mathbf{B}_0$. We discuss kinetic, density and
magnetic field fluctuations from the pencil beams, and identify physical
sources of non-Gaussian components to the PDF as single, strong shocks coupled
to fast magnetosonic compressions that form along $\mathbf{B}_0$. We discuss
the Gaussianisation of the $\mathcal{M} \gtrsim 4$ $s$-fields through the lens
of two phenomenologies: the self-similarity of the $s$-field and homogenisation
of the dynamical timescales between the over- and under-dense regions in the
compressible gas.
