With the rapid development of animal phenomics and deep phenotyping, we can get thousands of traditional but also molecular phenotypes per individual. However, there is still a lack of exploration regarding how to handle this huge amount of data in the context of animal breeding, presenting a challenge that we are likely to encounter more and more in the future. This study aimed to (1) explore the use of the Mega-scale linear mixed model (MegaLMM), a factor model-based approach, able to simultaneously estimate (co)variance components and genetic parameters in the context of thousands of milk traits, hereafter called thousand-trait (TT) models; (2) compare the phenotype values and genomic breeding values (u) predictions for focal traits (i.e., traits that are targeted for prediction, compared with secondary traits that are helping to evaluate), from single-trait (ST) and TT models, respectively; (3) propose a new approximate method of estimated genomic breeding values (U) prediction with TT models and MegaLMM. 3,421 milk mid-infrared (MIR) spectra wavepoints (called secondary traits) and 3 focal traits [average fat percent (Fat), average methane (CH4), and average somatic cell score (SCS)] collected on 3,302 first-parity Holstein cows were used. The 3,421 milk MIR wavepoints traits were composed of 311 wavepoints in 11 classes (months in lactation). Genotyping information of 564,439 SNP was available for all animals and was used to calculate the genomic relationship matrix. The MegaLMM was implemented in the framework of the Bayesian sparse factor model and solved through Gibbs sampling (Markov chain Monte Carlo). The heritabilities of the studied 3,421 milk MIR wavepoints gradually increased and then decreased in units of 311 wavepoints throughout the lactation. The genetic and phenotypic correlations between the first 311 wavepoints and the other 3,110 wavepoints were low. The accuracies of phenotype predictions from the ST model were lower than those from the TT model for Fat (0.51 vs. 0.93), CH4 (0.30 vs. 0.86), and SCS (0.14 vs. 0.33). The same trend was observed for the accuracies of u predictions: Fat (0.59 vs. 0.86), CH4 (0.47 vs. 0.78), and SCS (0.39 vs. 0.59). The average correlation between U predicted from the TT model and the new approximate method was 0.90. The new approximate method used for estimating U in MegaLMM will enhance the suitability of MegaLMM for applications in animal breeding. This study conducted an initial investigation into the application of thousands of traits in animal breeding and showed that the TT model is beneficial for the prediction of focal traits (phenotype and breeding values), especially for difficult-to-measure traits (e.g., CH4)., (The Authors. Published by Elsevier Inc. on behalf of the American Dairy Science Association®. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).)