Variational two-electron reduced-density-matrix (2RDM) theory is applied to computing energy spectra and properties of few-electron quantum dots. The model system is a two-dimensional electron gas with a central confinement potential. For each orbital angular momentum $J$, the energy and 2RDM are computed by the variational 2RDM method in which the energy is minimized as a functional of the 2RDM. In the minimization, which is performed by semidefinite programming, the 2RDM is constrained to represent a $N$-electron wave function with angular momentum $J$ by $N$- and $J$-representability conditions [D. A. Mazziotti, Phys. Rev. Lett. 93, 213001 (2004)]. Advantages of the variational 2RDM method include (i) lower bounds on the energies for all $J$ values, (ii) calculation of approximate 2RDMs in polynomial time without many-electron wave functions, (iii) exploitation of angular symmetry in the sparse block-diagonal structure of the 2RDM, (iv) accurate description of multireference correlation (entanglement) effects, and (v) direct calculation of one- and two-electron properties from the 2RDM. With the 2RDM we directly compute pair-correlation functions, radial charge densities, and average radial electron displacements in the quantum dot. Energies and properties are compared to those from solving the $N$-electron Schr\odinger equation by large-scale, exact diagonalization. It is found that the accuracy of the variational 2RDM approach is sensitive to the total orbital angular momentum and the symmetry of the final wave function. For quantum dots of high symmetry, the variational algorithm isolates a highly accurate solution that recovers the correlation energy within a few percent.