In this paper we focus on strong solutions of some heat-like problems with a non-local derivative in time induced by a Bernstein function and an elliptic operator given by the generator or the Fokker-Planck operator of a Pearson diffusion. Such kind of non-local equations naturally arise in the treatment of particle motion in heterogeneous media. In particular, we use spectral decomposition results for the usual Pearson diffusion to exploit explicit solutions of the aforementioned equations. Moreover, we provide stochastic representation of such solutions in terms of time-changed Pearson diffusions. Finally, we exploit some further properties of these processes, such as limit distributions and long/short-range dependence.