Direct Statistical Simulation (DSS) solves the equations of motion for the statistics of turbulent flows in place of the traditional route of accumulating statistics by Direct Numerical Simulation (DNS). That low-order statistics usually evolve slowly compared with instantaneous dynamics is one important advantage of DSS. Depending on the symmetry of the problem and the choice of averaging operation, however, DSS is usually more expensive computationally than DNS because even low order statistics typically have higher dimension than the underlying fields. Here we show that it is possible to go much further by using Proper Orthogonal Decomposition (POD) to address the "curse of dimensionality." We apply POD directly to DSS in the form of expansions in the equal-time cumulants to second order (CE2). We explore two averaging operations (zonal and ensemble) and test the approach on two idealized barotropic models on a rotating sphere (a jet that relaxes deterministically towards an unstable profile, and a stochastically-driven flow that spontaneously organizes into jets). Order-of-magnitude savings in computational cost are obtained in the reduced basis, potentially enabling access to parameter regimes beyond the reach of DNS.