We present a search for the production of neutral Higgs bosons in association with bottom quarks in $p ar p$ collisions at $sqrt s = 1.96$ TeV. The cross section for this process is enhanced in many extensions of the Standard Model (SM), such as in its Minimal Supersymmetric extension (MSSM) at large tan $eta$. The data, corresponding to a recorded integrated luminosity of 1 $fb^{-1}$, were collected with the DO detector at the Fermilab Tevatron Collider. In the absence of a signal a $95\%$ C.L. limit is set on the production cross section times branching ratio, and the results are also interpreted in the MSSM.