The nearly incompressible theory of magnetohydrodynamics (MHD) is formulated in the presence of a static large-scale inhomogeneous background. The theory is an inhomogeneous generalization of the homogeneous nearly incompressible MHD description of Zank & Matthaeus and a polytropic equation of state is assumed. The theory is primarily developed to describe solar wind turbulence where the assumption of a composition of two-dimensional (2D) and slab turbulence with the dominance of the 2D component has been used for some time. It was however unclear, if in the presence of a large-scale inhomogeneous background, the dominant component will also be mainly 2D and we consider three distinct MHD regimes for the plasma beta beta << 1, beta similar to 1, and beta >> 1. For regimes appropriate to the solar wind (beta << 1, beta similar to 1), compared to the homogeneous description of Zank & Matthaeus, the reduction of dimensionality for the leading-order description from three dimensional (3D) to 2D is only weak, with the parallel component of the velocity field proportional to the large-scale gradients in density and the magnetic field. Close to the Sun, however, where the large-scale magnetic field can be considered as purely radial, the collapse of dimensionality to 2D is complete. Leading-order density fluctuations are shown to be of the order of the sonic Mach number O(M) and evolve as a passive scalar mixed by the turbulent velocity field. It is emphasized that the usual "pseudosound" relation used to relate density and pressure fluctuations through the sound speed as delta rho = c(s)(2)delta p is not valid for the leading-order O(M) density fluctuations, and therefore in observational studies, the density fluctuations should not be analyzed through the pressure fluctuations. The pseudosound relation is valid only for higher order O(M-2) density fluctuations, and then only for short-length scales and fast timescales. The spectrum of the leading-order density fluctuations should be modeled as k(-5/3) in the inertial range, followed by a Bessel function solution K-nu(k), where for stationary turbulence nu = 1, in the viscous-convective and diffusion range. Other implications for solar wind turbulence with an emphasis on the evolution of density fluctuations are also discussed.