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 CD is computed as a function of the gas volume fraction φ and the Reynolds number Re = 2RUslip/ν for 0.03 φ 0.5 and 5 Re 2000. Here R, Uslip and ν 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 CD(Reeff)-relation of the suspension and the CD(Re)-relation of an individual bubble, where the effective Reynolds number Reeff = 2RUslip/νeff is based on the effective viscosity νeff which depends on the properties of the suspension. For Re ≲ 100, we find νeff ≈ ν/(1−0.6φ1/3), which is in qualitative agreement with previous proposed correlations for CD in bubble suspensions. For Re ≳ 100, on the other hand, we find νeff ≈ RUslipφ, which is explained by considering the turbulent kinetic energy levels in the liquid phase. Based on these findings, a correlation is constructed for CD(Re, φ). 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).