Variance components and breeding values for protein yield were estimated with REML without and with correction for heterogeneity of variances. Three different sire models were applied, which all accounted for genotype × environment (G × E) interaction. The first model included a sire × herd-year-season subclass (HYS) interaction. The second model divided all records in four different types of management groups, based on estimated HYS subclass effect. The third model, the reaction norm model, performed a random linear regression on the estimated HYS effect. For comparison, a standard model that did not take G × E interaction into account was also applied. Data consisted of 102,899 305-d first-lactation protein records of Holstein Friesians of 1000 of the largest Dutch dairy herds. All animals calved in 1997, 1998, or 1999. Estimated breeding values (EBV) for 2150 bulls with at least five daughters were calculated. The interaction model detected an interaction variance of 2.5% of the phenotypic variance. The EBV showed a correlation of 1.00 with those of the standard model without interaction. The model with the division in groups showed correlations between groups ranging from 0.73 to 0.86. The EBV showed correlations from 0.84 to 0.91 with the EBV of the standard model. The reaction norm model calculated EBV that had a correlation of 1.00 with the EBV of the standard model. The reaction norm model was not able to detect significant variance of the slope for the protein data corrected for heterogeneity of variances.