We extend the recently-developed theory of bulk orbital magnetization to finite electric fields, and use it to calculate the orbital magnetoelectric response of periodic insulators. Working in the independent-particle framework, we find that the finite-field orbital magnetization can be written as a sum of three gauge-invariant contributions, one of which has no counterpart at zero field. The extra contribution is collinear with and explicitly dependent on the electric field. The expression for the orbital magnetization is suitable for first-principles implementations, allowing to calculate the magnetoelectric response coefficients by numerical differentiation. Alternatively, perturbation-theory techniques may be used, and for that purpose we derive an expression directly for the linear magnetoelectric tensor by taking the first field-derivative analytically. Two types of terms are obtained. One, the `Chern-Simons' term, depends only on the unperturbed occupied orbitals and is purely isotropic. The other, `Kubo' term, involves the first-order change in the orbitals and gives isotropic as well as anisotropic contributions to the response. In ordinary magnetoelectric insulators both terms are generally present, while in strong Z2 topological insulators only the former type is allowed, and is quantized. In order to validate the theory we have calculated under periodic boundary conditions the linear magnetoelectric coupling for a 3-D tight-binding model of an ordinary magnetoelectric insulator, using both the finite-field and perturbation-theory expressions. The results are in excellent agreement with calculations on bounded samples.