A new four-field system of equations is derived from the compressible magnetohydrodynamic (MHD) equations for low Mach number turbulence in the solar wind and the interstellar medium, permeated by a spatially varying magnetic field. The plasma beta is assumed to be of order unity or less. It is shown that the full MHD equations can be reduced rigorously to a closed system for four fluctuating field variables: magnetic flux, vorticity, pressure, and parallel flow. Although the velocity perpendicular to the magnetic field is shown to obey a two-dimensional incompressibility condition (analogous to the Proudman-Taylor theorem in hydrodynamics), the three-dimensional dynamics exhibit the effects of compressibility. In the presence of spatial inhomogeneities, the four dynamical equations are coupled to each other, and pressure fluctuations enter the weakly compressible dynamics at leading order. If there are no spatial inhomogeneities and or if the plasma beta is low, the four-field equations reduce to the well-known equations of reduced magnetohydrodynamics (RMHD). For pressure-balanced structures, the four-field equations undergo a remarkable simplification which provides insight on the special nature of the fluctuations driven by these structures. The important role of spatial inhomogeneities is elucidated by 2.5-dimensional numerical simulations. In the presence of inhomogeneities, the saturated pressure and density fluctuations scale with the Mach number of the turbulence, and the system attains equipartition with respect to the kinetic, magnetic, and thermal energy of the fluctuations. The present work suggests that if heliospheric and interstellar turbulence exists in a plasma with large-scale, nonturbulent spatial gradients, one expects the pressure and density fluctuations to be of significantly larger magnitude than suggested in nearly incompressible models such as pseudosound.