We propose a model order reduction approach for non-intrusive surrogate modeling of parametric dynamical systems. The reduced model over the whole parameter space is built by combining surrogates in frequency only, built at few selected values of the parameters. This, in particular, requires matching the respective poles by solving an optimization problem. If the frequency surrogates are constructed by a suitable rational interpolation strategy, frequency and parameters can both be sampled in an adaptive fashion. This, in general, yields frequency surrogates with different numbers of poles, a situation addressed by our proposed algorithm. Moreover, we explain how our method can be applied even in high-dimensional settings, by employing locally-refined sparse grids in parameter space to weaken the curse of dimensionality. Numerical examples are used to showcase the effectiveness of the method, and to highlight some of its limitations in dealing with unbalanced pole matching, as well as with a large number of parameters.