Introduction
The number of secretory cells in the mammary gland is critically important to milk production. As a proxy for secretory cell numbers, the mammary size was associated with milk yield in goats, sheep and cattle, and the correlation coefficients between milk yield and the measures of mammary size (udder weight, udder volume and udder capacity) were moderate (Davis, Reference Davis2017). Moreover, udder conformation traits are often cited as being critically important indicators because mammary morphology is highly associated with the length of productive lif and mastitis susceptibility. Therefore, it is important to understand drivers of mammary gland morphology and their roles in the determination of productive performance in dairy livestock.
Genetic selection has increased udder size in Holstein cow and Friesian cows. Genome-wide association studies (GWAS) have facilitated identification of genes and variants associated with udder size. In Fleckvieh cattle, a total of 12 quantitative trait loci (QTL) and eight candidate genes were associated with tested variants (such as udder depth and fore udder length), and among the candidate genes, SP5 and ADAM12 may participate in regulation of mammary development (Pausch et al., Reference Pausch, Emmerling, Schwarzenbacher and Fries2016). Moreover, Flury (Reference Flury2014) reported a genome-wide association study for 13 udder morphology traits using imputed high-density genotypes in Brown Swiss sires. Five significant signals that affect different morphological features of mammary gland were detected, as well as one region located on chromosome six which was associated with clinical mastitis and milking speed. However, no genetic characteristics were reported for buffalo mammary gland morphology.
In this study, we performed genome-wide association study with 90 K Affymetrix Axiom Buffalo SNP Array to identify candidate genes and QTL that affect different morphological features of buffalo mammary gland. Five traits were used as response variables in 176 genotyped buffaloes. Our results are potentially beneficial for understanding the genetic regulation of buffalo mammary gland morphology.
Materials and methods
Animals and phenotypic records
The Ethical Animal Care and Use Committee of Federico II University of Naples (Italy) approved the experimental design and animal treatment.
A total of 176 Italian Mediterranean buffaloes born from 2005 to 2011 within at least one lactation were used for this study. The buffalo mammary gland morphology, including anterior teat length (ATL), posterior teat length (PTL), distance between anterior teats (AAD), distance between anterior and posterior teat (APD) and distance between posterior teats (PPD), were measured after milking, and every index were measured twice at different times of 1 d.
Genotyping and quality control
The DNAs were genotyped using the 90 K Axiom Buffalo SNP Array at Delta Genomics (Edmonton AB, Canada). Subsequently, quality control was performed using PLINK software version 1.9. Single nucleotide polymorphisms (SNPs) were excluded from the data set if they exhibited: (1) animal call rate <97%; (2) SNP call rate <95%; (3) minor allele frequency <0.05; (4) Hardy–Weinberg equilibrium (HWE) test P-value <10−6; and (5) chromosome and position were unknown.
Single trait genome wide association studies
Buffalo with both phenotypic data and genotype data were used to performed GWAS by using mixed model as following by both PLINK 1.07 and ASReml 2.0.
$$y_{ijklm}\,{\rm}= \,{\rm \mu}+{\rm parit}{\rm y}_i{\rm}+{\rm L}{\rm P}_j{\rm}+{\rm ag}{\rm e}_k{\rm}+{\rm SN}{\rm P}_l{\rm anima}{\rm l}_m{\rm}+ e_{ijklm}$$where y is the vector of phenotypes; μ is overall mean; parityi is fixed effect of the ith parity; LPj is the fixed effect of the jth lactation period (early, middle and late); agek is the kth covariate of age at calving; SNPl is fixed effect of the lth genotype of the SNP; animalm is mth random animal effect and e is the random residual error. Fixed effects were only fitted if they made a significant contribution to the overall variance component (P < 0.05). False discovery rate (FDR) was used to adjust P value for all available SNPs, and the significant threshold of SNP was proposed as FDR <0.1. The suggestive threshold was proposed as P < 10−4.
Haplotype analyses
The linkage disequilibrium (LD) and haplotype frequencies were obtained by Haploview 4.2 software. Haplotype blocks were constructed by using PHASE software, and the association between each haplotype block and mammary gland morphology was estimated in R. We used suggestive P value (P < 0.05) as threshold for considering candidate haplotype blocks affecting buffalo mammary gland morphology.
Results
Genome-wide association studies
Descriptive statistical results of the five indicators were presented in Supplementary File Table S1. After filtering low quality data by quality control, a total of 63 605 available SNPs and 176 buffaloes remained for subsequent GWAS. With single marker and single trait association analysis between buffalo mammary gland morphology and SNPs for each trait (Supplementary File Fig. S1), a total of 79 suggestive SNPs (P < 10−4) associated with at least one of five mammary gland morphology were identified (Supplementary File Table S2). After FDR adjustment, three SNPs (AX-85117983, AX-85094752, and AX-85117518) were significantly associated with buffalo ATL, PTL and APD, respectively. According to the map information, the three mutations were located on chromosome 1, chromosome 10, and chromosome 20, respectively.
Candidate genomic regions
To further identify QTL associated with buffalo udder traits, linkage disequilibrium (LD) analyses were performed to determine the relationship of SNPs within 0.5Mb window around three significant SNPs. A total of three haplotype blocks were recognized (Fig. 1). For each block, the association with buffalo mammary gland morphology was evaluated (Table 1). The block_2 was found to have an effect on buffalo PTL (P = 0.041), and block_3 was showed to be associated with buffalo ATL (P = 0.03) as well as PPD (P = 0.014). However, the block_1 were showed to be no significant effect on all buffalo mammary gland morphology in the present study.
Fig. 1. Haplotype blocks pattern for the significant SNPs based on the linkage disequilibrium within detected regions (0.5 Mb). Linkage disequilibrium block is presented in triangle box based on LD (squared correlation coefficient, r 2).
Table 1. The association between haplotypes and buffalo mammary gland morphology
ATL, anterior teat length; PTL, posterior teat length; AAD, distance between anterior teats; APD, distance between anterior and posterior teat; PPD, distance between posterior teats.
*P < 0.05: significant association after Bonferroni multiple test.
Candidate genes within identified genomic region
A total of 6 genes were discovered within the identified genomic region. Among the genes, GJD2, ACTC1, AQR, and ZNF770 located in the block_2 region, IPO11 and DIMT1 in the block_3 region. The description of these genes is presented in Supplementary File Table S3.
Discussion
In the present study, two significant blocks located on the 10th chromosomes and 20th chromosome respectively showed effects on buffalo udder traits. Interestingly, the areas where the two blocks were located were reported to be within or near the QTLs of mastitis (Talstein et al., Reference Talstein, Fontanesi, Dolezal, Scotti, Bagnato, Russo, Canavesi, Friedmann, Soller and Lipkin2010) and production life (Chamberlain et al., Reference Chamberlain, Hayes, Savin, Bolormaa, McPartlan, Bowman, Van der Jagt, MacEachern and Goddard2012). Furthermore, the block on chromosome 20 detected in this study is in a genomic region that controls milk yield traits in cattle (Arranz et al., Reference Arranz, Coppieters, Berzi, Cambisano, Grisart, Karim, Marcq, Moreau, Mezer and Riquet2015). However, we did not detect any SNP associated with buffalo udder traits on chromosome 5, which is a chromosome that is recognized to include many QTLs related to udder traits in dairy cattle and beef breeds.
Endocrine factors play a crucial role in the regulation of mammary growth (Rowson et al., Reference Rowson, Daniels, Ellis and Hovey2012). In the present study, a total of 6 genes were located on the significant region, and one endocrine-related gene was included. GJD2 (Gap Junction Protein Delta 2), located on chromosomes 10 from 30.255 to 30.257 Mb, encodes a member of the connexin protein family, which may be involved in the regulation of insulin secretion (Valentina et al., Reference Valentina, Celine, Pierri, Samuel, Jacques-Antoine, João, Valeriya, Leif, Rico and Fabrizio2016). Insulin was considered to play a role in mammogenesis and lactation, since the concentration of insulin is related to the expression of leptin, which was detected in adipocytes, alveolar epithelial, and mammary epithelial cells (Çatlı et al., Reference Çatlı, Olgaç Dündar and Dündar2014; Palin et al., Reference Palin, Farmer and Duarte2017). Therefore, we speculate that GJD2 gene may affect mammogenesis in buffalo by regulating expression of insulin and leptin. Having said this, it should be pointed out that the major mammogenic endocrine signals comprise steroid hormones and polypeptide hormones of pituitary and placental origin.
ACTC1 (Cardiac alpha-actin) located on the block_2 region was identified to be associated with buffalo mammary gland morphology. ACTC1 belongs to the actin family and plays a crucial role in the muscle formation of filaments. The polymorphisms of ACTC1 were associated with hypertrophic cardiomyopathy (Despond and Dawson, Reference Despond and Dawson2018). Two studies of selection signatures in dairy and beef cattle revealed many interesting candidate genes affecting milk production, body size and muscle formation, and both of them included ACTC1 (Qanbari et al., Reference Qanbari, Gianola, Hayes, Schenkel, Miller, Moore, Thaller and Simianer2011; Zhao et al., Reference Zhao, McParland, Kearney, Du and Berry2015). Moreover, ACTC1 expressed in the retina and implicated as a candidate gene that influence the development of myopia (Solouki et al., Reference Solouki, Verhoeven, van Duijn, Verkerk, Ikram, Hysi, Despriet, van Koolwijk, Ho and Ramdas2010). Interestingly, a gene associated with progressive vision loss was reported to be associated with udder structure in cattle breeds (Michot et al., Reference Michot, Chahory, Marete, Grohs, Dagios, Donzel, Aboukadiri, Deloche, Allais-Bonnet, Chambrial, Barbey, Genestout, Boussaha, Danchin-Burge, Fritz, Boichard and Capitan2016), although the biological relevance of this observation is not clear.
The association between udder traits and milk production performance has been broadly corroborated, and many GWAS reported candidate genes may influence mammary development (Flury, Reference Flury2014; Pausch et al., Reference Pausch, Emmerling, Schwarzenbacher and Fries2016). However, most of these genes have a limited role in milk productive traits through effects on mammary growth, and many genes associated with milk productive traits were not found to be associated with udder traits, such as DGAT1 (Li et al., Reference Li, Liu, Li, Zhang, Hua, Salzano, Campanile, Gasparrini, Liang and Yang2018), PRL (Li et al., Reference Li, Liang, Li, Du, Hua, Salzano, Campanile, Gasparrini and Yang2017), and ABCG2 (Cohen-Zinder et al., Reference Cohen-Zinder, Seroussi, Larkin, Loor, Everts-van, Lee, Drackley, Band, Hernandez, Shani, Lewin, Weller and Ron2005). In the present study, none of candidate genes except ACTC1 have obvious effects on milk production.
In conclusion, this study revealed candidate genes and potential genomic regions for mammary gland morphology in buffalo. To our knowledge, this is the first report about genetic basis of buffalo mammary gland morphology by genome-wide association study. These findings will provide a better understanding of the genetic regulation of mammary growth. However, development and lactation of mammary gland are complex biological process, and more data are needed to answer the question whether genetic improvement in mammary gland morphology will be beneficial for the productive performance.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0022029919000967.
Acknowledgements
This work was supported by the Earmarked Fund for Modern Agri-industry Technology Research System (CARS-37-04B). We thanks for all field staff and participants in our study.
