Combined effects of CXCL8 (IL-8) and CXCR2 (IL-8R) gene polymorphisms on deregressed MACE EBV indexes of milk-related traits in Simmental bulls

Abstract CXCL8 (also known as IL-8) is a member of the CXC subfamily of chemokines that binds two of the seven transmembrane G-protein-coupled receptors (GPCRs), CXCR1 and CXCR2, to mediate and regulate leucocyte accumulation and activation at sites of inflammation. They are known to play a critical role in both disease susceptibility and infection outcome. The aim of this study was to investigate the entire sequences of CXCL8 and CXCR2 genes in thirty-one Simmental sires to evaluate the effects of genomic variants on the indexes of the bulls for milk, fat and protein yields, and for somatic cell score (SCS). Five new single nucleotide polymorphisms (SNPs) were found in CXCR2 gene. The analysis of association indicated that one SNP in CXCL8 and two in CXCR2 influenced the considered traits. To evaluate the existence of functional haplotypic effects, combinations among the three genomic variants (SNP 1 in CXCL8, SNP 6 and SNP 7 in CXCR2) were investigated. Four different haplotypic alleles were identified in the experimental population, one of which at a high frequency (61%). Bulls with Hap 4 (G-C-G at SNP 1, SNP 6, and SNP 7 respectively) had more favourable indexes for SCS (P < 0.05). These results suggest that the SNPs in CXCL8 and CXCR2 may be potential genetic markers to improve udder health in the Simmental breed.

CXCL8 (IL8) is a proinflammatory chemokine belonging to the C-X-C family, whose functions extend to innate and adaptive immune cell lineages, thereby indicating a critical role in both disease susceptibility and infection outcome (Mukaida, 2000). CXCL8 is produced by several cell types including epithelial cells, and signals through the ligation with CXCR1 and CXCR2, two members of the seven-transmembrane G protein-coupled receptor family, which has CXCR2 as its primary functional receptor (Liu et al., 2016). Binding of CXCL8 to its receptors on the neutrophil surface induces neutrophil activation, stimulates chemotaxis and increases phagocytosis and killing ability (Mitchell et al., 2003).
In Bos taurus, the CXCL8 gene is located on chromosome 6 with a transcript reported in GenBank Accession Number: NM_173925 and a coded protein composed of 101 amino acids (NP_776350). Several CXCL8 SNPs have been identified (Heaton et al., 2001) and associated with milk fat yield and udder depth (Leyva-Baca et al., 2007;Chen et al., 2011). Furthermore, Meade et al. (2012) identified several (29) polymorphic sites across a 2.1 kb upstream promoter region of CXCL8 and the sequence analysis identified two distinct promoter haplotypes (IL8-h1 and IL8-h2). Significant inter-breed differences in haplotype frequencies were found in Holstein-Friesian, Norwegian Red, Jersey and Italian Simmental breeds (Meade et al., 2012;Stojkovic et al., 2016;De Matteis et al., 2021).
The C-X-C motif chemokine receptor 2 (CXCR2) gene is located on chromosome 2 and two isoforms X1 (GenBank Accession Number: XM_024978406) and X2 (XM_024978410) are transcripted. The coded proteins are made of 400 (XP_024834174) and 387 (XP_024834178) amino acids, respectively. Single nucleotide polymorphisms within the CXCR2 gene have been identified and evaluated for their potential association with disease in humans (Kato et al., 2000;Renzoni et al., 2000;Yang, et al., 2006;Javor et al., 2012) and in cattle (Rambeaud and Pighetti, 2005). In bovine, Youngerman et al. (2004) identified five single nucleotide polymorphisms in the CXCR2 gene. Rambeaud and Pighetti (2005) showed that the SNP16 (+777 G→C) (Table 1) results in amino acid substitution and Holstein cows carrying the CC genotype had an increased incidence of subclinical mastitis compared to cows that expressed the CG or GG genotype. Moreover, Beecher et al. (2010) showed that the G allele of SNP16 tended to associate with decreased somatic cell score (SCS) throughout the entire lactation, as well as with increased fat yield. This paper documents the presence of a substantial number of polymorphisms in CXCL8 and CXCR2 genes, some of which combined in functional haplotypes are associated with SCS, an indicator of udder health.

Animals and data
The entire CXCL8 and CXCR2 gene sequences of 31 Simmental sires were investigated to evaluate the effects of genomic variants on their deregressed multiple across country evaluation (MACE) estimated breeding values (EBV) for milk, fat and protein yields, and for SCS. The experimental design procedure was detailed in our previous study involving the same groups of animals (Napolitano et al., 2021). All data were provided by the Italian Simmental Breeders Association (ANAPRI). Sequences from thirty-one Simmental sires from different origins (20 Austrian, 6 Swiss, 3 German and 2 Italian), born between 1981 and 2006, with an average of 6316 (95-42 657) daughters were analysed. Indexes concerning yields of milk, milk fat, milk protein and SCS refer to a 3-year period (2013)(2014)(2015).

Genetic variants in the CXCL8 and CXCR2 genes
In order to investigate the presence of polymorphisms in the coding and regulatory regions of CXCL8 and CXCR2 genes, we used the sequences of the genomic regions of chromosome 6 (GenBank, Accession Number: NC_037333 for CXCL8) and chromosome 2 (GenBank, Accession Number: NC_037329 for CXCR2) in Dominette's Bos taurus. These sequences range between the nucleotides g.88810001…88815000 (CXCL8) and g.106184001… 106195000 (CXCR2), containing both the coding and regulatory regions of these genes. These chromosomal regions were blasted with the whole genomic sequencing of each of the 31 bulls to highlight any polymorphisms. The same sequences were checked through the Ensembl archive (https://www.ensembl.org/info/website/index.html) to verify if the identified genomic variants had already been reported. Genome-wide sequencing of each bull was uploaded to Galaxy server at https://usegalaxy.eu (Version 2.3.4.3) and analysed as described by Napolitano et al. (2021). All the polymorphisms identified in the CXCL8 and CXCR2 genes were individually tested to evaluate their influence on the productive traits examined. The ones that produced a significant effect were then evaluated together using a functional haplotype.

Statistical analysis
On each SNP site, the χ 2 test for deviation from Hardy-Weinberg equilibrium, along with expected and observed heterozygosity, polymorphism information content (PIC) and linkage disequilibrium (LD) were calculated using the algorithms provided by SAS software 9.4 (ALLELE procedure).
The associations between each SNP and milk related traits were analysed using the general linear model (GLM procedure of SAS software 9.4): each genotype was independently modelled as a fixed factor in Y ijn = μ + G j + e ijn , where Y ijn is the phenotype for trait i of animal n carrying the genotype j, μ the overall population mean; G j the fixed effect of the genotype j, with j = 1, 2, 3 depending if the animal n is homozygous for one allele, heterozygous, or homozygous for the other allele; and e ijn the random error.
In order to evaluate associations between functional haplotypes and milk related traits, each haplotype was independently modelled as a fixed factor in Y ijkn = μ + H jk + e ijkn , where Y ijkn is the phenotype of trait i for animal n carrying the haplotype j in k copies, μ the overall population mean; H jk the fixed effect of haplotype j, with k = 1, 2, 3 depending if the animal n carries none, one or two copies of the haplotype j; and e ijkn the random error.
The statistical significance of all traits and least-square means were determined by Tukey's test available in the GLM procedure (LSMEANS/ADJUST = TUKEY).

Diversity analysis of CXCL8 and CXCR2 SNPs
Twenty-eight SNPs were identified in the described genomic sequences (Table 1) twenty-five of them causing sequence variations in coding or regulatory regions of the analysed genes. Among those positioned in coding regions, SNPs 9, 13, 14 and 16 (exon 2 of the CXCR2), are missense mutations (GenBank, Accession Number: XP_024834174) while the others are synonymous (8, 10-12, 15 and 17-21 of exon 3 of the CXCR2). Four out of the new five SNPs were found in CXCR2 exon 3 (SNP 8,9,10 and 19).
Almost all SNPs identified on the two examined genes were informative markers (online Supplementary Table S1), with a PIC value over 0.20. Genotypes detected at the markers of IL8 were in Hardy-Weinberg equilibrium, contrary to what was found on its receptor where most of the markers were not (χ 2 = 6.05-21.19; P < 0.01 to 0.0001).
Correlations among SNPs are reported in online Supplementary Table S2. In the CXCL8 gene, SNP 1 was 100% in linkage disequilibrium (LD = 1) with SNPs 4 and 5 and in 95% LD with SNPs 2 and 3. Fifteen SNPs of the CXCR2 gene with LD < 0.9: SNPs 6, 7, 9 (novel), 11, SNPs 13-18, 19 (novel), 20, 21, 23, and 27 were processed for the estimation of their associations with phenotypes (deregressed MACE EBV indexes of milk, fat and protein yields, and SCS) as a function of their genotype. Only SNP 1 was chosen from the CXCL8 gene for the association analysis based on its location in the regulatory region of the gene, as previously done for association tests in dairy cattle by Meade et al. (2012). Table 2 reports the association of the genotypes at the 16 SNP loci with the indexes on milk, fat and proteins yields as well as SCS. Association analysis indicated that only six out of the sixteen chosen markers affected the considered traits. The SNP1 of CXCL8 gene exhibited a significant effect (P < 0.05) on all the evaluated indexes. In the CXCR2 gene, SNP6 showed a significant effect on milk yield, protein yield and SCS, SNP7 influenced milk yield and SCS indexes, SNPs 9 and 17 influenced milk yield, while SNP13 influenced fat yield, all at P < 0.05.

Association analysis of functional haplotypes on milk indexes
As reported in Table 2, only 3 out of the 16 analysed SNPs (1, 6 and 7) were determined to have significant effects on SCS. These 3 variants were investigated together using a functional haplotype.
Functional haplotype make-up and frequencies are listed in Table 3. Four different haplotypic alleles with a frequency ranging from 8 to 61% were identified. Association analyses between functional haplotypes in CXCL8 and CXCR2 genes with milk related indexes are shown in Table 4. Haplotype Hap1 (ACG) influenced fat and protein yield as well as SCS score. Hap2 (GAA) influenced milk and milk protein yield, and SCS score, but with a significant difference (P < 0.05) only when the haplotype occurred in homozygosis. The haplotype Hap4 (GCG), the most frequent in our population (61%), improved the SCS index (P < 0.05) leaving the quantitative aspects of the milk unchanged.

Discussion
We investigated in a cohort of thirty-one Simmental sires, the possible effects of genomic variants in CXCL8 and CXCR2 on the deregressed MACE EBV indexes for quality and quantity of milk, as well as SCS score, an indicator of udder health status. The identified variants in both genes were located in regulatory regions (5 ′ and 3 ′ UTR) and within intronic and exonic sequences. Many of these polymorphisms, as already reported by other authors (Youngerman et al., 2004;Meade et al., 2012), were in perfect linkage disequilibrium in the surveyed 31 genomes.
The most interesting SNP in CXCL8 was SNP1 (g.88810697 G > A), previously identified by Meade et al. (2012) and used in association tests in dairy cattle. SNP1 is located 5 bp upstream from the NFκB (C-rel) binding site and occurs within the TFBS for C/EBP and NFAT transcription factors. Such proximity of SNP1 to the NF-κB binding site suggests a potentially important regulatory role in bovine CXCL8 gene expression. The 'A' allele of SNP1 introduces two predicted TFBS for Oct-1 that would be abrogated by the 'G' allele (Meade et al., 2012). The Oct-1 transcriptional repressor can repress CXCL8 expression (Sibbet et al., 1995;Bhat et al., 1996;Zhang et al., 1999) by displacing the C/EBP transcription enhancer from the CXCL8 gene promoter (Wu et al., 1997). Based on the human model (dela Paz et al., 2007), removal of the Oct-1 repressor binding site would be expected to upregulate IL-8 production in cattle possessing 'G' allele. Meade et al. (2012) found two different haplotypes (IL8-h1 and IL8-h2) of the IL8 gene. The 'A' allele is carried by the IL8-h1 haplotype, whereas the 'G' allele by the IL8-h2 one. These authors evaluated the implications of both haplotypes for the bovine immune response and demonstrated that cows carrying the IL8-h2 showed a higher IL8 protein expression at 12 h post in vivo LPS stimulation compared to IL8-h1 haplotype. The practical relevance of this haplotype was the detection of a genetic association between IL8-h2 and somatic cell counta marker of mastitis (Stojkovic et al., 2017). Furthermore, divergent cattle populations with different selection pressures for health and production traits showed different frequencies for the 'G' and 'A' alleles (Meade et al., 2012;Stojkovic et al., 2016;De Matteis et al., 2021).
In the present study, we observed 68% of Simmental bulls carrying homozygous allele 'G' at SNP1 and only one sire with homozygous allele 'A'. Moreover, the association analysis showed that allele 'A' highlights significant negative effect on quantitative and qualitative traits as well as on SCS score. This result agrees with previous association studies showing a positive effect of 'G' allele to increase SCS (Stojkovic et al., 2017). Other previous studies have already reported SNPs within the bovine CXCR2 gene (Rambeaud and Pighetti, 2005) associated with somatic cell score (Leyva-Baca et al., 2008;Goertz et al., 2009) and mastitis resistance (Youngerman et al., 2004) as well as with other production, health and reproductive traits (Galvao et al., 2011). These results support the importance of the IL-8 axis in disease resistance and productive lifespan in dairy cattle.
In this study, five new SNPs were identified in the CXCR2 gene and one of these, SNP9, resulted in amino acid 89 histidinetyrosine substitution. One of the previously identified variations in CXCR2 gene, the SNP16 (g.106192177 G > C ) resulted in amino acid 245 (285aa on XP_024978406 protein sequence) glutamine-histidine replacement (Rambeaud and Pighetti, 2005). This polymorphism affects receptor activation, because the region is important for G-protein coupling and activation (Damaj et al., 1996). Subsequent studies showed that Holstein cows expressing the CC genotype at position +777 (SNP16) had an increased incidence of subclinical mastitis compared to Holstein cows that expressed the CG or the GG genotype (Youngerman et al., 2004). In addition, they exhibited impaired neutrophil migration and adhesion molecule upregulation compared to cows with GG genotype. It is noteworthy that no CC genotype was found in our sample of Simmental bulls and no significant difference was detected on the considered milk traits and SCS between GG and GC genotypes.
Three out of 16 analysed SNPs (1 in the CXCL8 gene, 6 and 7 in the CXCR2 gene) produced significant effects on SCS and at least on one milk related index. Therefore, to evaluate the existence of haplotypic effects, SNP1, SNP6 and SNP7 were used to  (2019) proposed a new method based on functional haplotype which considers both the main and epistatic effects among SNPs, to overcome some constraints of the GWAS in which only consecutive SNPs were evaluated and, therefore, only additive and dominance effects were considered. Compared with single SNP analysis, the combined genotype analysis provides more information on gene interactions. This supports the notion that a relatively large number of variables at functionally relevant loci exert their influence on complex trait variation primarily via epistatic interactions, rather than through conventional additive and dominance effects (Jarvis and Cheverud, 2011). The functional haplotypic approach allowed us to gain a systemic view of the probable role of the CXCL8 and the CXCR2 in the determinism of milk traits. In our study, the combination effects of CXCL8-SNP1, CXCR2-SNP6 and CXCR2-SNP7 significantly affected quantitative and qualitative milk traits and SCS score. Among the combined haplotypes, the Hap4 was the most frequent (61%) and bulls with homozygous G-C-G genotype combination showed the highest, most favourable, SCS score, corresponding to the lowest somatic cell count (SCC). These results suggest that the SNPs in CXCL8 and CXCR2 may be potential genetic markers for SCS not only for the Simmental breed but also for other dairy breeds. For example, it should be possible to detect in advance the frequency of haplotype Hap4 (G-C-G) in the population, if the selection goals include milk quality and mastitis resistance. Unfortunately, homozygous and heterozygous Hap4 bulls tended to have the lowest milk, fat and protein productions, even if the differences were not significant. Even if results need to be confirmed at a larger scale, haplotype effects showed the same pattern: a higher SCS (i.e. lower SCC) results in a lower milk yield. The trend was consistent throughout the four haplotypes with some variability due to few observations in some classes.
In conclusion, this study, within the limits of our sample size, implements knowledge on CXCL8 and CXCR2 genomic variability in Simmental breed highlighting five new variants in CXCR2. Furthermore, the multiple locus analysis revealed that combined effects of CXCL8 and CXCR2 are likely to affect SCS and other milk-related deregressed MACE EBV indexes. The significant association of CXCL8-CXCR2 functional haplotypes with SCS score supports the hypothesis that this genotype combination may have important functional implications for the expression of IL8 and ultimately on bovine immune response, particularly in mammary epithelial cells. However, studies on a larger scale are needed to verify if these genes could be used in actual breeding programmes to increase dairy cow resilience.