Abstract We investigate finite-temperature observables in three-dimensional large N critical vector models taking into account the effects suppressed by $$ \frac{1}{N} $$ 1 N . Such subleading contributions are captured by the fluctuations of the Hubbard-Stratonovich auxiliary field which need to be handled with care due to a subtle divergence structure which we clarify. The examples we consider include the scalar O(N) model, the Gross-Neveu model, the Nambu-Jona-Lasinio model and the massless Chern-Simons Quantum Electrodynamics. We present explicit results for the free energy density to the subleading order in $$ \frac{1}{N} $$ 1 N , which captures the thermal one-point function of the stress-energy tensor to this order. We also include the dependence on a chemical potential. We determine the Wilson coefficient in the thermal effective action that is sensitive to global symmetry for the first time directly in interacting CFTs, which produces a symmetry-resolved asymptotic density of states. We further provide a formula from diagrammatics for the one-point functions of general single-trace higher-spin currents. We observe that in most cases considered, these subleading effects lift the apparent degeneracies between observables in different models at infinite N, while in special cases the discrepancies only start to appear at the next-to-subleading order.