We investigate the problem of deriving highest density regions (HDRs) from multivariate data samples. We are interested in estimating minimum volume sets that contain a given probability. In the case of unknown distribution probabilities f, the problem involves their estimation, which may be challenging in multidimensional settings. Motivated by the ubiquitous role of copula modelling in modern statistics, we explore their use in the context of HDR estimation. Rather than directly estimating the multivariate f, we propose to estimate the marginals and their dependence structure, i.e., the copula structure, separately. We evaluate this new method, using both a parametric and a nonparametric approach, in a number of synthetic experiments and considering a real dataset.