In multiple-traits animal models, variance components are estimated using REstricted Maximum Likelihood (REML) method. Best Linear Unbiased Prediction (BLUP) of breeding values is then performed and used for selection. Such an approach is based on the assumption of normality of each bivariate phenotypes even if, in practice, this hypothesis is not always realistic. Indeed, in a multivariate setting, even if the marginal distributions of each phenotypes are normally distributed, the multi-normality assumption may be violated due to the non-normality for the cross-sectional dependence structure, that is to say the copula of the multivariate distribution fails to be Gaussian. We assessed by simulation, the impact of a non-Gaussian copula for the residual term of the mixed model on the REML estimations of the genetic parameters (heritability, correlation) and on the BLUP of breeding values. Large populations of 18510 animals over 6 generations were simulated and bivariate phenotypes were generated considering different copulas (symmetric, asymmetric, with tail dependence or not) for the dependence structure between the error components. Results obtained highlighted the robustness of REML estimation of genetic and environmental variance-covariance components and prediction of breeding values in a multivariate mixed model context when the Gaussian hypothesis is not verified.