In recent years, neural quantum states have emerged as a powerful variational approach, achieving state-of-the-art accuracy when representing the ground-state wave function of a great variety of quantum many-body systems, including spin lattices, interacting fermions or continuous-variable systems. However, accurate neural representations of the ground state of interacting bosons on a lattice have remained elusive. We introduce a neural backflow Jastrow Ansatz, in which occupation factors are dressed with translationally equivariant many-body features generated by a deep neural network. We show that this neural quantum state is able to faithfully represent the ground state of the 2D Bose-Hubbard Hamiltonian across all values of the interaction strength. We scale our simulations to lattices of dimension up to 20×20 while achieving the best variational energies reported for this model. This enables us to investigate the scaling of the entanglement entropy across the superfluid-to-Mott quantum phase transition, a quantity hard to extract with non-variational approaches.