Markovian population models are suitable abstractions to describe well-mixed interacting particle systems in situation where stochastic fluctuations are significant due to the involvement of low copy particles. In molecular biology, measurements on the single-cell level attest to this stochasticity and one is tempted to interpret such measurements across an isogenic cell population as different sample paths of one and the same Markov model. Over recent years evidence built up against this interpretation due to the presence of cell-to-cell variability stemming from factors other than intrinsic fluctuations. To account for this extrinsic variability, Markovian models in random environments need to be considered and a key emerging question is how to perform inference for such models. We model extrinsic variability by a random parametrization of all propensity functions. To detect which of those propensities have significant variability, we lay out a sparse learning procedure captured by a hierarchical Bayesian model whose evidence function is iteratively maximized using a variational Bayesian expectation-maximization algorithm.