The study was carried out on three Polish Black and White Holstein-Friesian dairy herds (A–C). Herds were considered to be among the best in Poland in terms of MY. A total of 782 dairy cows were used in the study. Cows were fed a total mixed ration and had access to feed ad libitum.

Milk recording data were collected in the period from 2016 to 2019 according to A4 method (the interval between two successive recordings of two daily milkings ranging between 28 and 33 days) accredited by the International Committee for Animal Recording (ICAR) [14], by the Milk Analysis Laboratory of the Polish Federation of Cattle Breeders and Dairy Farmers (PFCB&DF) which is certified by ISO 17025 and certified by ICAR. All data were stored and provided by SYMLEK IT system (by PFCB&DF). Genetic data used in this study were collected in the process of routine estimation of breeding value (EBV) that is performed in Poland. Poland as a member of the Euro Genomics cooperative uses customised EuroGenomics arrays. In the present study the following versions were used: Eurogenomics v3_POL; Eurogenomics v4_POL; Eurogenomics v5_POL; Eurogenomics v6_POL; Eurogenomics v8b_POL; Eurogenomics MD_POL with Infinium HD Illumina protocol. The different panels had very little overlap and due to limited number of genotypes available the imputation was not successful, thus only 6,765 genome wide SNPs with call rate >0.95 across all animals and minor allele frequency >0.001 were used.

Biopsy samples (ear punch) were collected in the period from 2015 to 2019. Eurogenomics arrays’ probe list and sequences are confidential and could not be published in the present paper. For the purpose of the present study authors have chosen SNPs that were present on all array versions. Farmers have authorised the authors to access a portion of the SNP data. Genotyping data provided by PFCB&DF and farmers were evaluated and 10 polymorphisms within PAEP (in exons 2, 3, 4, and 5) were identified (Figure 1, Table 1).

Allele and genotype frequencies for all SNPs were calculated and data was subjected to the Hardy-Weinberg equilibrium test. The national milk recording scheme (SYMLEK system) provided test day records for somatic cell counts (SCC, ths/mL) in the milk of all studied cows. SCC data were log transformed to the natural logarithm for somatic cell count (LnSCC).

A dataset containing 9,284 records on milk performance was analysed. The following parameters were recorded during the research: MY (kg/d), percentage of fat (FP, %), protein (PP, %), dry matter (DMP, %) and lactose (LP, %), urea content (UC, mg/L) as well as natural LnSCC (ln). Moreover, a dataset containing information on PAEP genotypes was added into the analysis.

In order to check the validity of the use of haplotypes for GEBV estimation, two models were used: with and without haplotypes in pedigree based and single-step calculations. The 7-trait model (repeatability) was used with fixed effects of herd, year-season of the sampling, year-season of calving, sample number (related to days in milk), sample number squared and age of cow at sampling; with the last three effects fitted as covariates. For this purpose, free Fortran 90 software was used [15]. A total of 191 cows born in 2017 was selected for the validation group. The method of partial vs whole data was used based on the research carried out by Legarra and Reverter [16] who showed that the correlation of EBVs in validation animals between the whole and partial data set reflects ratio of the accuracies of the two models. The model was run with pedigree as the relationship matrix and as a single step model. Analysed traits were statistically characterised. The impact of respective PAEP composite genotypes on tested parameters was assessed by one-way analysis of variance followed by the Bonferroni post-test. p-Values less than 0.05 were considered as statistically significant (SAS software package) [17].

In the present study 10 genetic variants (SNPs) were investigated and allele and genotype frequencies were calculated (Table 2). It was noted that seven genotypes were monomorphic, while three (PAEP_3982, PAEP_5174, PAEP_5263) showed genetic diversity. CT genotype was the most frequent amongst PAEP_3982 genotypes (50%), with CC and TT being less frequent (30% and 20%, respectively). In our study, Hardy-Weinberg equilibrium test was performed showing observed and expected heterozygosities at the same levels (0.30), chi² at the level of 0.208. p Value (0.649) was greater than α = 0.05, which indicates that, regarding PAEP_3982 SNP, the studied population was consistent with Hardy-Weinberg equilibrium. Similarly, 49% of all tested cows were heterozygous in terms of PAEP_5174, while 32% had TT and 19% CC genotypes. Population was in Hardy-Weinberg equilibrium as chi2 (0.0145) and p value (0.904) did not allow this hypothesis to be disregarded. Results showed that 49% of all individuals were heterozygous for PAEP_5263, the Hardy-Weinberg equilibrium test revealed that the population was also consistent with Hardy-Weinberg equilibrium (p = 0.904). In the present study, seven out of ten SNPs were monomorphic, which may suggest that due to extensive breeding and selection, genetic diversity for these SNPs was reduced. Other authors also have investigated PAEP variants, however most focused mainly on variants A and B, which are characterised by a number of different SNPs [6,12,18], among others three SNPs that were detected to be polymorphic in the present study. Comparing the data showing nucleotide presence in each variant, it can be noted that variant A is characterised by the presence of thymine in PAEP_3982, cytosine in PAEP_5174 and thymine in PAEP_5263, while variant B has cytosine in the position 3,982, thymine in 5,174 and cytosine in 5,263. Hristov et al [12] reported that the frequencies of A and B alleles differ according to cattle breed: for Shorthorn Rhodopean cattle breed AA and BB frequencies were low (0.063) and heterozygote frequency was very high (0.88), while in Bulgarian Rhodopean cattle the frequency of AA genotype was higher (0.40) and AB lower (0.58). Also, Barbosa et al [9] investigated polymorphisms within the PAEP gene revealing two alleles and three genotypes. They reported that within the population of Girolando cows, most common were heterozygotes (0.48) and noted that genotype frequencies were in Hardy-Weinberg equilibrium (chi² = 1.04; p<0.01).

Moreover, out of the three PAEP variants (PAEP_3982, PAEP_5174, PAEP_5263) that showed polymorphism, 5 composite genotypes were formed (Table 3): P1 (PAEP_3982CC +PAEP_5174TT+PAEP_5263CC, giving the composite genotype CCTTCC), P2 (CTCTCT), P3 (CTTTCC), P4 (TTCCTT), P5 (TTCTCT). The most frequent composite genotype in the studied population was P2, with 48% animals having that composite genotype, while P3 and P5 were least frequent (these composite genotypes were found only in 2% of the tested population). Further studies on a larger population may further help with analysing the frequencies of created composite genotypes.

The relative increase in correlation calculated based on partial vs whole data by including haplotypes in the GEBV calculation is presented in Table 4. In 5 out of 7 analysed traits, some improvement in accuracy of the EBVs by incorporating haplotypes was found with an average of 15.8% for pedigree based EBVs and 9.5% for single step GEBVs increase in accuracy across the traits. Overall, the accuracy was low because the population structure was not favourable and the genetic relationship between training and validation was very low, also the amount of available data used in the training was small. With the relative gains for MY, UC, and SCC in the range of 30% for pedigree-based analysis, genotyping three SNPs necessary to determine haplotypes may be beneficial for breeding programmes with limited resources that can not afford the use of microarrays. However, validation on a larger dataset and in the target population is recommended.

Significant differences (p<0.05) were found between different PAEP composite genotypes regarding tested traits. Differences between composite genotype groups (P1 to P5) in relation to different parameters are presented in Table 5. Animals with P5 composite genotype were characterised by the highest (p<0.05) values for PP and LnSCC and the lowest for LP and UC (p<0.05). P4 composite genotype was linked with increased DMP and UC, while P3 was associated with increased LP as well as decreased PP and LnSCC. While no statistical differences were noted while studying tested herds as one group, differences in MY between composite genotypes were noted in Herd A and B. These results are in accordance with our previous reports; for instance, Sitkowska et al [19] who investigated polymorphism in the PAEP gene reported that PAEP genotypes have an impact on MY and milk composition. They noted that animals with BB genotype were characterised by a higher FP and PP in milk, while allele A was associated with MY. Similarly, Dokso et al [13] noticed that MY of animals with AA genotype was higher by 9.58% than for BB genotype, while cows homozygous for B allele had a higher FP (by 4.28%) than AA animals. Wageh Zaglool et al [4] linked AA genotype to a higher MY and PP (by 16.81% compared to BB animals) and BB with higher FP in milk (by 30.95% compared to AA animals). These authors did not report significant differences in LP between different genotypes.

Monitoring UC in milk of dairy cattle may help farmers to detect health issues that may result in reduced fertility of their cows. Concentration of milk urea provides information on the efficiency of protein synthesis, as well as on the balance between crude protein and energy. As urea dissolves in milk, milk urea concentration is strongly correlated with blood urea nitrogen level. A high level of UC suggests an inefficient protein conversion. On the other hand, low UC may cause health and fertility problems [20]. Literature suggests that the optimal UC in cow’s milk is between 150 and 300 mg/L [20,21]. In the present study UC did not exceed 300 mg/L (on the average UC was 245.01 to 270.1 mg/L) (Table 5). Lower values were reported by Czajkowska et al [21] (206.4 mg/L) and Rzewuska and Strabel [20] (230.1 mg/L).

The SCC is an important parameter that has a great economic impact since higher levels of somatic cells are associated with mastitis [22]. In the present study, a significant association between PAEP composite genotypes and LnSCC was recorded in the whole population as well as in herd A and C, indicating that animals with P5 composite genotype were characterised by the highest and P3 by the lowest LnSCC (Table 5). In herd B, the same tendency was observed, however, the differences were not statistically significant.

It is worth noticing that selected herds for this study were among the best in Poland and where farmers have incorporated genomic selection into their breeding programmes for many years. Their work may contribute to the fact that MY obtained by tested cows was higher than many mean values for the average dairy cattle herd. For instance, while authors of the present study reported MY to be 37.69 to 38.40 kg/d, the average MY in 2018 in Poland was 27.21 kg/d [23]. Mulder et al [24] noted MY to be even lower: 24.31 kg/d. Soyudal et al [6], who investigated the influence of respective PAEP alleles on milk performance of Holstein Friesian cows in Turkey, reported MY in 305-day lactation to be 28.90 kg/d for animals with genotype AA, 27.86 kg/d kg for heterozygotes and 26.36 kg/d for BB. These levels were considerably lower than those reported in the present study. The differences between MYs in different herds may result from many factors, including different management. Herds selected in the present study are considered as one of the best in Poland in terms of MY, and the owners have put an emphasis on a careful breeding management, that was also supported by genotyping data.

Microarray genotyping in dairy herds is used mostly for EBV, the data on specific SNPs is often not further investigated. Confirming the impact of selected PAEP SNPs that are included in the microarrays and their haplotypes may allow farmers to obtain information on genotypes related to only these SNPs that are relevant for their production systems. It is worth mentioning that many farmers decide not to genotype their herds with microarrays due to high costs. EBV performed based on pedigrees combined with genetic markers evaluation is more accurate than EBV that does not include information on genetic markers. Therefore confirming the impact of PAEP SNPs and their haplotypes on milk performance may allow farmers to obtain relevant genetic information in more cost-efficient way as genotyping three SNPs would be cheaper than genotyping with the use of microarray. Fitting the haplotypes as fixed effects in addition to the presence of PAEP SNPs in the construction of the H matrix also provided some gains in accuracy which suggests that the effect of these haplotypes is larger than assumed by a single step GBLUP model.