  We study the ground state properties of spin-1 bosons in a two-dimensional
optical lattice, by applying a variational Monte Carlo method to the S=1
Bose-Hubbard model on a square lattice at unit filling. A doublon-holon binding
factor introduced in the trial state provides a noticeable improvement in the
variational energy over the conventional Gutzwiller wave function and allows us
to deal effectively with the inter-site correlations of particle densities and
spins. We systematically show how spin-dependent interactions modify the
superfluid-Mott insulator transitions in the S=1 Bose-Hubbard model due to the
interplay between the density and spin fluctuations of bosons. Furthermore,
regarding the magnetic phases in the Mott region, the calculated spin structure
factor elucidates the emergence of nematic and ferromagnetic spin orders for
antiferromagnetic ($U_2>0$) and ferromagnetic ($U_2<0$) couplings,
respectively.
