We study nonlinear convection in a rapidly rotating sphere with internal heating for values of the Prandtl number relevant for liquid metals (Pr ∈ [10 −2 , 10 −1 ]). We use a numerical model based on the quasi-geostrophic approximation, in which variations of the axial vorticity along the rotation axis are neglected, whereas the temperature field is fully three-dimensional. We identify two separate branches of convection close to onset: (i) a well-known weak branch for Ekman numbers greater than 10 −6 , which is continuous at the onset (supercritical bifurcation) and consists of thermal Rossby waves, and (ii) a novel strong branch at lower Ekman numbers, which is discontinuous at the onset. The strong branch becomes subcritical for Ekman numbers of the order of 10 −8 . On the strong branch, the Reynolds number of the flow is greater than 10 3 , and a strong zonal flow with multiple jets develops, even close to the nonlinear onset of convection. We find that the subcriticality is amplified by decreasing the Prandtl number. The two branches can co-exist for intermediate Ekman numbers, leading to hysteresis (Ek = 10 −6 , Pr = 10 −2 ). Non-linear oscillations are observed near the onset of convection for Ek = 10 −7 and Pr = 10 −1 .