We use numerical simulations to address locomotion at zero Reynolds number in viscoelastic (Giesekus) fluids. The swimmers are assumed to be spherical, to self-propel using tangential surface deformation, and the computations are implemented using a finite element method. The emphasis of the study is on the change of the swimming kinematics, energetics, and flow disturbance from Newtonian to viscoelastic, and on the distinction between pusher and puller swimmers. In all cases, the viscoelastic swimming speed is below the Newtonian one, with a minimum obtained for intermediate values of the Weissenberg number, We. An analysis of the flow field places the origin of this swimming degradation in non-Newtonian elongational stresses. The power required for swimming is also systematically below the Newtonian power, and always a decreasing function of We. A detail energetic balance of the swimming problem points at the polymeric part of the stress as the primary We-decreasing energetic contribution, while the contributions of the work done by the swimmer from the solvent remain essentially We-independent. In addition, we observe negative values of the polymeric power density in some flow regions, indicating positive elastic work by the polymers on the fluid. The hydrodynamic efficiency, defined as the ratio of the useful to total rate of work, is always above the Newtonian case, with a maximum relative value obtained at intermediate Weissenberg numbers. Finally, the presence of polymeric stresses leads to an increase of the rate of decay of the flow velocity in the fluid, and a decrease of the magnitude of the stresslet governing the magnitude of the effective bulk stress in the fluid.