In this paper we introduce a class of time-changed processes obtained by composing a Pearson diffusion with the inverse of a subordinator. Such time-changed processes provide stochastic representations of solutions of some Cauchy problems with a non-local derivative in time induced by a suitable Bernstein function. In particular we show the existence of strong solutions for these equations via spectral methods and then we use the spectral decomposition results to deduce first order stationary and limit distributions of the introduced processes. Finally, we show that the first order stationary processes are not second order stationary, neither in wide-sense, and they exhibit short or long range dependence (in some sense) depending on the Bernstein function.