Lattice Boltzmann and immersed boundary methods are used to conduct direct numerical simulations of suspensions of massless, spherical gas bubbles driven by buoyancy in a three-dimensional periodic domain. The drag coefficient C(D) is computed as a function of the gas volume fraction phi and the Reynolds number R(e) = 2RU(slip)/nu for 0.03 (sic) phi (sic) 0.5 and 5 (sic) Re (sic) 2000. Here R, U(slip) and nu denote the bubble radius, the slip velocity between the liquid and the gas phases and the kinematic viscosity of the liquid phase, respectively. The results are rationalized by assuming a similarity between the C(D)(Re(eff))-relation of the suspension and the C(D)(Re)-relation of an individual bubble, where the effective Reynolds number Re(eff) = 2RU(slip)/nu(eff) is based on the effective viscosity nu(eff) which depends on the properties of the suspension. For Re less than or similar to 100, we find nu(eff) approximate to nu/(1-0.6 phi(1/3)), which is in qualitative agreement with previous proposed correlations for C(D) in bubble suspensions. For Re greater than or similar to 100, on the other hand, we find nu(eff) approximate to RU(slip)phi, which is explained by considering the turbulent kinetic energy levels in the liquid phase. Based on these findings, a correlation is constructed for C(D)(Re, phi). A modification of the drag correlation is proposed to account for effects of bubble deformation, by the inclusion of a correction factor based on the theory of Moore (J. Fluid Mech., vol. 23, 1995, p. 749).