In this paper the numerical approximation of stochastic differential equations satisfying a global monotonicity condition is studied. The strong rate of convergence with respect to the mean square norm is determined to be $1/2$ for the two-step BDF-Maruyama scheme and for the backward Euler-Maruyama method. In particular, this is the first paper which proves a strong convergence rate for a multi-step method applied to equations with possibly superlinearly growing drift and diffusion coefficient functions. We also present numerical experiments for the $3/2$-volatility model from finance and a two dimensional problem related to Galerkin approximation of SPDE, which verify our results in practice and indicate that the BDF2-Maruyama method offers advantages over Euler-type methods if the stochastic differential equation is stiff or driven by a noise with small intensity.