For m number of bosons, carrying spin (S=1) degree of freedom, in Ω number of single particle orbitals, each triply degenerate, we introduce and analyze embedded Gaussian orthogonal ensemble of random matrices generated by random two-body interactions that are spin (S) scalar [BEGOE(2)-S1]. The embedding algebra is U (3) ⊃ G ⊃ G1 ⊗ SO(3) with SO(3) generating spin S. A method for constructing the ensembles in fixed-(m, S) space has been developed. Numerical calculations show that the form of the fixed-(m, S) density of states is close to Gaussian and level fluctuations follow GOE. Propagation formulas for the fixed-(m, S) space energy centroids and spectral variances are derived for a general one plus two-body Hamiltonian preserving spin. In addition to these, we also introduce two different pairing symmetry algebras in the space defined by BEGOE(2)-S1 and the structure of ground states is studied for each paring symmetry.