We study convexity properties of the average integral operators naturally associated with divergence-form second-order subelliptic operators ℒ with nonnegative characteristic form. When ℒ is the classical Laplace operator, these average operators are the usual average integrals over Euclidean spheres. In our subelliptic setting, the average operators are (weighted) integrals over the level sets ∂ Ω r ( x ) = { y : Γ ( x , y ) = 1 ∕ r } of the fundamental solution Γ(x,y) of ℒ. We shall obtain characterizations of the ℒ-subharmonic functions u (that is, the weak solutions to −ℒu≤0) in terms of the convexity (w.r.t. a power of r) of the average of u over ∂Ωr(x), as a function of the radius r. Solid average operators will be considered as well. Our main tools are representation formulae of the (weak) derivatives of the average operators w.r.t. the radius. As applications, we shall obtain Poisson–Jensen and Bocher type results for ℒ.