Tolerance to heat stress is one of the most relevant traits in defining the ability of animals to cope with environmental challenges. Phenotypes for measuring heat tolerance could be obtained by combining production and climate variables. Principal Component Analysis (PCA) carried out on test day records previously corrected for systematic effects except for THI is able to derive 2 principal components (PC) that summarize the individual patterns of corrected production data across different THI levels. However PCA does not allow for missing data. To overcome this problem, in this work 35,992 TD records of for milk yield, fat percentage and somatic cell score (SCS) measured on Italian Holstein cows sired by 3,697 bulls were analyzed by a mixed linear model that included the fixed effects of calving month nested within year, age at calving, herd nested within year, and of the interaction between days in milk and parity, and the random effect of the sire. Sire (co)variance was modeled by 11 × 11 unstructured matrix, where 11 were the classes of THI into which TD records were grouped. PCA was then carried out on the estimated matrices. Finally, PC scores for each trait were calculated by multiplying the eigenvector matrix for the matrix of standardized data. As expected, the first principal component (PC1) was a measure of the level at which the pattern of corrected TD records across THI level is located, whereas the second (PC2) summarized the shape of this curve. Heritabilities of the PC1 scores were 0.39, 0.44, and 0.37 for milk yield, fat percentage, and SCS respectively. Values for PC2 were 0.07, 0.05, and 0.06 respectively. Genetic correlations between PC1 and PC2 were lower than 0.20 for milk yield and fat percentage, 0.39 for SCS. The combined use of mixed linear model and PCA allowed the derivation of 2 phenotypes able to summarize describe the overall level and the shape of patterns of milk production traits across different levels of THI index.