A variety of swimming microorganisms, called ciliates, exploit the bending of a large number of small and densely packed organelles, termed cilia, in order to propel themselves in a viscous fluid. We consider a spherical envelope model for such ciliary locomotion where the dynamics of the individual cilia are replaced by that of a continuous overlaying surface allowed to deform tangentially to itself. Employing a variational approach, we determine numerically the time-periodic deformation of such surface which leads to low-Reynolds locomotion with minimum rate of energy dissipation (maximum efficiency). Employing both Lagrangian and Eulerian points of views, we show that in the optimal swimming stroke, individual cilia display weak asymmetric beating, but that a significant symmetry-breaking occurs at the organism level, with the whole surface deforming in a wavelike fashion reminiscent of metachronal waves of biological cilia. This wave motion is analyzed using a formal modal decomposition, is found to occur in the same direction as the swimming direction, and is interpreted as due to a spatial distribution of phase differences in the kinematics of individual cilia. Using additional constrained optimizations, as well as a constructed analytical ansatz, we derive a complete optimization diagram where all swimming efficiencies, swimming speeds, and amplitudes of surface deformation can be reached, with the mathematically optimal swimmer, of efficiency one-half, being a singular limit. Biologically, our work suggests therefore that metachronal waves may allow cilia to propel cells forward while reducing the energy dissipated in the surrounding fluid.