We investigate quarkonium mass spectra in external constant magnetic fields by using QCD sum rules. We first discuss a general framework of QCD sum rules necessary for properly extracting meson spectra from current correlators computed in the presence of strong magnetic fields, that is, a consistent treatment of mixing effects caused in the mesonic degrees of freedom. We then implement operator product expansions for pseudoscalar and vector heavy-quark current correlators by taking into account external constant magnetic fields as operators, and obtain mass shifts of the lowestlying bound states ηc and J/ψ in the static limit with their vanishing spatial momenta. Comparing results from QCD sum rules with those from hadronic effective theories, we find that the dominant origin of mass shifts comes from a mixing between ηc and J/ψ with a longitudinal spin polarization, accompanied by other subdominant effects such as mixing with higher excited states and continua.