In this paper we extend the local iterative Lie-Schwinger block-diagonalization method - introduced in [DFPR3] for quantum lattice systems with bounded interactions in arbitrary dimension- to systems with unbounded interactions, i.e., systems of bosons. We study Hamiltonians that can be written as the sum of a gapped operator consisting of a sum of on-site terms and a perturbation given by relatively bounded (but unbounded) interaction potentials of short range multiplied by a real coupling constant t. For sufficiently small values of |t| independent of the size of the lattice, we prove that the spectral gap above the ground-state energy of such Hamiltonians remains strictly positive. As in [DFPR3], we iteratively construct a sequence of local block-diagonalization steps based on unitary conjugations of the original Hamiltonian and inspired by the Lie-Schwinger procedure. To control the ranges and supports of the effective potentials generated in the course of our block-diagonalization steps, we use methods introduced in [DFPR3] for Hamiltonians with bounded interactions potentials. However, due to the unboundedness of the interaction potentails, weighted operator norms must be introduced, and some of the steps of the inductive proof by which we control the weighted norms of the effective potentials require special care to cope with matrix elements of unbounded operators. We stress that no large-field problems appear in our construction. In this respect our operator methods turn out to be an efficient tool to separate the low-energy spectral region of the Hamiltonian from other spectral regions, where the unbounded nature of the interaction potentials would become manifest.