An approximation for treating multiple quantum nuclei within the nuclear-electronic orbital (NEO) framework for molecular systems is presented. In the approximation to NEO-Hartree-Fock, the nuclear wave function is represented by a Hartree product rather than a Slater determinant, corresponding to the neglect of the nuclear exchange interactions. In the approximation to NEO-density functional theory, the nuclear exchange-correlation functional is chosen to be the diagonal nuclear exchange interaction terms, thereby eliminating the nuclear self-interaction terms. To further enhance the simplicity and computational efficiency, the nuclear molecular orbitals or Kohn-Sham orbitals are expanded in terms of localized nuclear basis sets. These approximations are valid because of the inherent localization of the nuclear orbitals and the numerical insignificance of the nuclear exchange interactions in molecular systems. Moreover, these approximations lead to substantial computational savings due to the reduction in both the number of integrals that must be calculated and the size of the matrices that must be diagonalized. These nuclear Hartree product approximation (HPA) methods scale linearly with the number of quantum protons and are highly parallelizable. Applications to a water hexamer, glycine dimer, and 32-water cluster, where all hydrogen nuclei are treated quantum mechanically, illustrate the accuracy and computational efficiency of the nuclear HPA methods. These strategies will facilitate the implementation of explicitly correlated NEO methods for molecular systems with multiple quantum protons.