In multiple-trait animal models, variance components are frequently estimated using Restricted Maximum Likelihood method (REML). Such an approach assumes the multivariate normality for the phenotype even if, in practice, this hypothesis is not always realistic. We assessed, using simulation, the impact of a non-Gaussian distribution for the residual term of the mixed model, on the REML estimations. The non-Gaussian distributions were simulated using a copula-based approach. Large populations over 8 generations were simulated using random selection for the 3 first generations and using a truncation selection for the following. Results obtained highlighted the robustness of REML when random selection is performed. On the contrary with a truncation selection process, we observed significant differences with the true parameters, particularly with asymmetric bivariate distributions on the residual part.