A finite element code based on the level-set method is used to perform direct numerical simulations (DNS) of the transient and steady-state motion of bubbles rising in a viscoelastic liquid modelled by the Oldroyd-B constitutive equation. The role of the governing dimensionless parameters, the capillary number (Ca), the Deborah number (De) and the polymer concentration parameter c, in both the rising speed and the deformation of the bubbles is studied. Simulations show that there exists a critical bubble volume at which there is a sharp increase in the terminal velocity with increasing bubble volume, similar to the behaviour observed in experiments, and that the shape of both the bubble and its wake structure changes fundamentally at that critical volume value. The bubbles with volumes smaller than the critical volume are prolate shaped while those with volumes larger than the critical volume have cusp-like trailing ends. In the latter situation, we show that there is a net force in the upward direction because the surface tension no longer integrates to zero. In addition, the structure of the wake of a bubble with a volume smaller than the critical volume is similar to that of a bubble rising in a Newtonian fluid, whereas the wake structure of a bubble with a volume larger than the critical value is strikingly different. Specifically, in addition to the vortex ring located at the equator of the bubble similar to the one present for a Newtonian fluid, a vortex ring is also present in the wake of a larger bubble, with a circulation of opposite sign, thus corresponding to the formation of a negative wake. This not only coincides with the appearance of a cusp-like trailing end of the rising bubble but also propels the bubble, the direction of the fluid velocity behind the bubble being in the opposite direction to that of the bubble. These DNS results are in agreement with experiments.