On harmonic function spaces, we define shift operators using zonal harmonics and partial derivatives, and develop their basic properties. These operators turn out to be multiplications by the coordinate variables followed by projections on harmonic subspaces. This duality gives rise to a new identity for Gegenbauer polynomials. We introduce large families of reproducing kernel Hilbert spaces of harmonic functions on the unit ball of $\mathbb{R}^n$ and investigate the action of the shift operators on them. We prove a von Neumann inequality with harmonic polynomials for a subfamily of commuting row contractions that we call harmonic type. The right side of this inequality contains the norms of shift operators on a special harmonic Hilbert function space that we denote G ̆. Thus this space is a natural harmonic counterpart of the Drury-Arveson space. This is joint work with Daniel Alpay of Chapman University, Orange, CA.