<p>In quantitative genetics and breeding, phenotypes are measured to improve traits of interest. These phenotypes combine genetic and environmental components and are therefore modeled using mixed models. As genetic values are unobservable, they are treated as latent variables with a covariance structure defined by the pedigree. Multiple observed phenotypes are typically assumed to follow a joint Gaussian distribution. Estimation is then performed via restricted maximum likelihood (REML) or Bayesian methods under Gaussian assumptions. However, even if each phenotype’s component appears Gaussian, the joint distribution may exhibit lower tail dependence despite Gaussian marginals, due to complex dependencies between random variables. This can lead to biased estimates of variance components. We propose an extension of the standard genetic and environmental mixed model by incorporating copula functions to flexibly capture the joint distribution of phenotypes beyond Gaussian assumptions. We develop a stochastic gradient descent algorithm combined with a Monte Carlo Markov Chain step to estimate variance components and predict genetic values. The method is evaluated through simulations and applied to a real pig breeding dataset, where normality assumptions for the joint phenotype distribution seem inappropriate. Our results show that the copula-based approach provides more robust estimates compared to standard REML under non-Gaussian conditions.</p>