We introduce the Nuclear Electronic All-Particle Density Matrix Renormalization Group (NEAP-DMRG) method for solving the time-independent Schrodinger equation simultaneously for electrons and other quantum species. In contrast to already existing multicomponent approaches, in this work we construct from the outset a multi-reference trial wave function with stochastically optimized non-orthogonal Gaussian orbitals. By iterative refining of the Gaussians positions and widths, we obtain a compact multi-reference expansion for the multicomponent wave function. We extend the DMRG algorithm to multicomponent wave functions to take into account inter- and intra-species correlation effects. The efficient parametrization of the total wave function as a matrix product state allows NEAP-DMRG to accurately approximate full configuration interaction energies of molecular systems with more than three nuclei and twelve particles in total, which is currently a major challenge for other multicomponent approaches. We present NEAP-DMRG results for two few-body systems, i.e., H$_2$ and H$_3^+$, and one larger system, namely BH$_3$