Numerical studies on the weak Galerkin spectral element method for eigenvalue problems of second order elliptic equations are carried out in this paper. In analogy to weak Galerkin finite element methods, the approximation space in a weak Galerkin spectral element method contains independently interior components on subdomains, and boundary components at subdomain interfaces are supplemented to ensure information interchange between subdomains. The weak Galerkin quadrilateral spectral element method with arbitrarily convex quadrilateral meshes is the focus of this paper. The interior and boundary components of the weak approximation function on each subdomain are constructed by orthogonal polynomials on the reference square and its edges via the bilinear mapping, respectively. While approximation spaces for weak gradients on each subdomain are established by orthogonal polynomials on the reference square via the Piola transform. In such circumstances, a weak Galerkin quadrilateral spectral element approximation scheme together with its implementation algorithm is then proposed for eigenvalue problems of second order elliptic equations, and the wellposedness of the approximation scheme is analyzed in detail after a systematic investigation of nullities of discrete weak gradients. Abundant numerical experiments are performed to elaborately analyze the accuracy and convergence of the weak Galerkin quadrilateral spectral element method, especially the impact of the diverse match of polynomial degrees in the approximation function space and the discrete weak gradient space upon the accuracy and convergence. Numerical results show that the p-version weak Galerkin quadrilateral spectral element method inherits the exponential orders of convergence of spectral methods, while the h-version weak Galerkin quadrilateral spectral element method not only possesses full-order convergence of h-version methods but also achieves a super convergence by properly matching polynomial degrees in the approximation spaces.