Neural responses are variable: even under identical experimental conditions, single neuron and population responses typically differ from trial to trial and across time. Recent work demonstrated that this variability can be modulated by specific stimulus features, most prominently quenched by stimulus onset. However, it is unknown how ubiquitous this effect is across stimulus features, and indeed whether it exists in brain areas to which classical neural coding approaches cannot be applied due to the lack of trials with repeatable experimental stimuli. To address these questions, we develop a universal probabilistic spike count model that builds on sparse Gaussian processes and can model arbitrary spike count distributions (SCDs) with flexible dependence on observed covariates. Without requiring repeatable trials, our method can flexibly capture covariate-dependent joint SCDs, and thus provide a characterisation of neural variability and its modulation by a set of covariates. We apply the model to recordings from four different cortical and subcortical areas, including non-sensory neural populations: V1 (both mouse and monkey), head direction cells in the anterodorsal thalamic nucleus (ADn), hippocampal place cells, and entorhinal grid cells. We find that variability across brain areas defies a simple parametric relationship with mean spike count as assumed in standard models, and its modulation by a diverse set of external covariates can be comparably strong to that of the mean firing rate. These results demonstrate that the modulation of variability is a widespread phenomenon across the brain, with a diverse set of covariates contributing to it.