Let X be a multi-type continuous-state branching process with immigration on state space R+d. Denote by gt, t≥0, the law of X(t). We provide sufficient conditions under which gt has, for each t>0, a density with respect to the Lebesgue measure. Such density has, by construction, some Besov regularity. Our approach is based on a discrete integration by parts formula combined with a precise estimate on the error of the one-step Euler approximations of the process. As an auxiliary result, we also provide a criterion for the existence of densities of solutions to a general stochastic equation driven by Brownian motions and Poisson random measures, whose coefficients are Hölder continuous and might be unbounded.