Implications
The genomic regions subjected to selection tend to generate runs of homozygosity (ROH) islands or hotspots. The aim of this work was to identify the differences between breeds and use the location of ROH islands to identify genes potentially involved in economically important traits. We identified several genes within ROH involved in a wide variety of biological processes, such as milk yield and composition, reproduction, immune response, resistance/susceptibility to infectious and diseases. These results showed that scanning the genome for ROH could be an alternative strategy to detect genomic regions and genes related with important economically traits.
Introduction
The development of single nucleotide polymorphism (SNP) arrays to scan the genome allow us to distinguish non-autozygotic segments that are identical by state from autozygotic and identical by descent segments (Peripolli et al., Reference Peripolli, Munari, Silva, Lima, Irgang and Baldi2016). A potential alternative method, called ROH, has been used in livestock for the identification of homozygous genomic regions (Purfield et al., Reference Purfield, Berry, McParland and Bradley2012; Ferenčaković et al., Reference Ferenčaković, Hamzić, Gredler, Solberg, Klemetsdal, Curik and Sölkner2013a). Runs of homozygosity are contiguous lengths of homozygous genotypes that are present in an individual because the parents transmitted identical by descending haplotypes to their offspring (Gibson et al., Reference Gibson, Morton and Collins2006). Runs of homozygosity has been widely used as predictors of whole-genome inbreeding levels (Zhang et al., Reference Zhang, Calus, Guldbrandtsen, Lund and Sahana2015a; Mastrangelo et al., Reference Mastrangelo, Tolone, Di Gerlando, Fontanesi, Sardina and Portolano2016). Moreover, ROH have been used in livestock genomic studies, confirming the correlation between shared ROH and genomic regions putatively under selection (Kim et al., Reference Kim, Cole, Huson, Wiggans, Van Tassell, Crooker, Liu, Da and Sonstegard2013; Gaspa et al., Reference Gaspa, Marras, Sorbolini, Ajmone-Marsan, Williams, Valentini, Dimauro and Macciotta2014; Metzger et al., Reference Metzger, Karwath, Tonda, Beltran, Águeda, Gut, Gut and Distl2015; Szmatola et al., Reference Szmatola, Gurgul, Ropka-Molik, Jasielczuck, Zabek and Bugno-Poniewierska2016; Kukučková et al., Reference Kukučková, Moravčíková, Ferenčaković, Simčič, Mészáros, Sölkner, Trakovická, Kadlečík, Curik and Kasarda2017; Purfield et al., Reference Purfield, McParland, Wall and Berry2017). In fact, the genomic regions subjected to selection frequently show signatures, such as reduced nucleotide diversity, and tend to generate ROH islands or hotspots, which have high levels of homozygosity around a selected locus compared with the rest of the genome (Szmatola et al., Reference Szmatola, Gurgul, Ropka-Molik, Jasielczuck, Zabek and Bugno-Poniewierska2016; Purfield et al., Reference Purfield, McParland, Wall and Berry2017). Runs of homozygosity islands are not randomly distributed across the genome and are shared among individuals within a breed (Zhang et al., Reference Zhang, Guldbrandtsen, Bosse, Lund and Sahana2015b).
A large number of cattle breeds are defined by marked phenotypic differences and, therefore, constitute valuable models to study genome evolution in response to processes such as selection and domestication. Thus, in livestock species, ROH may contribute to the detection of genomic regions that could explain phenotypic differences among breeds that affect traits of economic importance. We previously described ROH structures in three local cattle breeds (Reggiana, Cinisara and Modicana) and in Holstein cattle (Mastrangelo et al., Reference Mastrangelo, Tolone, Di Gerlando, Fontanesi, Sardina and Portolano2016). The aim of this work was to further study the distribution of ROH islands across the genome of these four cattle breeds, which may provide insights into the mechanisms underlying their genomic differences. In addition, it aimed to characterize ROH islands and identify enriched genes that could potentially explain the effects of these homozygous regions on economically important traits.
Material and methods
Samples, genotyping and data filtering
A total of 407 animals (Cinisara=71, Modicana=72, Reggiana=168 and Italian Holstein=96) were used for the analyses. All of the individuals were genotyped using the Illumina BovineSNP50 v2 BeadChip assay (Illumina Inc., San Diego, CA, USA). Single nucleotide polymorphisms were filtered to exclude loci assigned to unmapped contigs, and only those SNPs located on autosomes were considered. Quality control included call frequency ⩾0.95, minor allele frequency (MAF) ⩾0.01, and Hardy–Weinberg Equilibrium with a P>0.001. SNPs that did not satisfy these quality criteria were excluded. Single nucleotide polymorphisms were mapped using the Bos taurus UMD 3.1.1 genome assembly.
Genetic relationship between individuals
The genetic relationship among individuals was estimated by principal components analysis (PCA) of genetic distances. This analysis was based on the identity by state (IBS) matrices of genetic distances between individuals. Principal components analysis of the genetic distance (D) matrix was performed using the multidimensional scaling option in PLINK v.1.07 (Purcell et al., Reference Purcell, Neale, Todd-Brown, Thomas, Ferreira, Bender, Maller, Sklar, de Bakker, Daly and Sham2007). The graphical representation was depicted using the statistical R software (http://www.R-project.org/).
Runs of homozygosity detection
Runs of homozygosity were estimated, for each individual, using a sliding window approach of 50 SNPs in PLINK v.1.07 (Purcell et al., Reference Purcell, Neale, Todd-Brown, Thomas, Ferreira, Bender, Maller, Sklar, de Bakker, Daly and Sham2007). The minimum length that constituted the ROH was set to 4 Mb. The density of the SNP panel used to generate data for ROH identification is an important factor that strongly affects autozygosity estimates. The 50K panel overestimates the number of small segments (Purfield et al., Reference Purfield, Berry, McParland and Bradley2012; Ferenčaković et al., Reference Ferenčaković, Solkner and Curik2013b). The following criteria were used to define the ROH: (i) one missing SNP was allowed in the ROH and up to one possible heterozygous genotype, (ii) the minimum number of consecutive SNPs that constituted a ROH was set to 30, (iii) minimum density of 1 SNP every 100 kb, and (iv) maximum gap between consecutive SNPs of 1 Mb.
Identification of genomic regions and genes within runs of homozygosity
To identify the genomic regions of high homozygosity, the amount of times that each SNP appeared in the ROH was considered and normalized by dividing it by the number of animals included in the analysis. These values were plotted against the position of the SNP along the chromosome. The genomic regions were defined according to Szmatola et al. (Reference Szmatola, Gurgul, Ropka-Molik, Jasielczuck, Zabek and Bugno-Poniewierska2016). Adjacent SNPs having a proportion of ROH occurrences over the adopted threshold formed ROH islands. Mean linkage disequilibrium (LD) was estimated using HAPLOVIEW v. 4.2 (Barrett et al., Reference Barrett, Fry, Maller and Daly2005) for all pairwise combinations of SNPs within each ROH island. Genomic coordinates for all identified ROH islands were also used for the annotation of genes that were fully or partially contained within each selected region using the UCSC Genome Browser (http://genome.ucsc.edu/). The genes were further analyzed with the Panther Classification System (Mi et al., Reference Mi, Muruganujan and Thomas2013) to identify significant (P⩽0.05) gene ontology (GO) terms. Finally, to investigate the biological function of each annotated gene contained in ROH islands, an accurate literature search was also conducted.
Results
A PCA was used to visualize and explore the genetic relationships among breeds. The PCA (Figure 1) showed that breeds formed non-overlapping clusters and were clearly separated populations. After data quality and genetic relationships analyses, no outliers were detected.
Genetic relationship defined with multidimensional scaling analysis for the four cattle breeds.
A total of 44 875 SNPs in Cinisara, 42 687 SNPs in Modicana, 35 270 SNPs in Reggiana, and 41 569 SNPs in Italian Holstein cattle breeds were retained after quality control for ROH detection. The top 1% of SNPs observed in the ROH was selected, and adjacent SNPs over this threshold were merged into genomic regions corresponding to ROH islands (Szmatola et al., Reference Szmatola, Gurgul, Ropka-Molik, Jasielczuck, Zabek and Bugno-Poniewierska2016). In ROH islands detected here, each SNP showed a percentage of occurrence >10% (Figure 2). This approach resulted in the identification of 11 ROH islands in Cinisara and Italian Holstein, and eight in Modicana and Reggiana (Table 1). Two overlapping ROH islands were observed between breed pairs. Modicana and Reggiana breeds showed a common genomic region on Bos taurus autosome (BTA) 6 (6:38 689 886 to 39 346 170 bp) and Cinisara and Italian Holstein breeds on BTA10 (10:56 464 919 to 56 792 715 bp). The genomic distribution of ROH islands was clearly non-uniform among breeds and across autosomes (Table 1). The longest ROH island was observed in Italian Holstein on BTA10 (12.42 Mb), while the shortest one was observed in Reggiana on BTA3 (0.03 Mb). BTA6 in Modicana breed had the ROH with the highest peak (Figure 2) which consisted of 38 SNPs with an occurrence in ROH >45% and a length of 2.05 Mb.
Genome-wide frequency of single nucleotide polymorphisms (SNPs) occurrence into runs of homozygosity (ROH) for each cattle breed. The black lines indicated the adopted threshold which defines the autozygosity islands, different per each breed (top 1% of the observations).
List of genomic regions of extended homozygosity (ROH islands) identified in each cattle breed
The mean r 2 value, a standard descriptive LD parameter, was estimated for all pairwise combinations of SNPs within each ROH island (Supplementary Material Table S1). In Cinisara breed, the majority of SNPs within ROH islands showed low level of LD (<0.080), and r 2 ranged from 0.024 to 0.290. The other breeds showed intermediate levels of LD within ROH islands (from 0.006 to 0.280). The highest LD level was found in the ROH island on BTA6 in Reggiana breed (0.387).
Within all of the ROH islands here reported, we identified from 126 to 347 genes (347 Italian Holstein, 250 Modicana, 190 Cinisara and 126 Reggiana). A list of genes found in the ROH islands of each breed underwent a GO enrichment analysis. Multiple categories were statistically significant (P⩽0.05). The genes within ROH islands encompass a wide spectrum of molecular function, biological process, and cellular components. A PANTHER gene list analysis revealed a high percentage of genes involved in catalytic activity (GO:0003824), cellular processes (GO:0009987), cell part (GO:0044464), metabolic processes (GO:0008152), binding (GO:0005488) as well as biological regulations (GO:0065007) and response to stimulus (GO:0050896) in all of the ROH islands of the analyzed breeds (Table 2). Supplementary Material Table S2 provides the chromosome position, number of SNPs and number of genes per genomic region, gene symbol and full name for all of the annotated genes in each breed.
Gene ontology (GO) terms enriched (P<0.05) based on runs of homozygosity islands and number of involved genes (n) for each cattle breed
Discussion
We analyzed animals from four Italian cattle breeds with different inbreeding background and selection histories. Mastrangelo et al. (Reference Mastrangelo, Tolone, Di Gerlando, Fontanesi, Sardina and Portolano2016), in a previous study on evaluation on ROH in these breeds, reported the highest value of inbreeding (F) based on ROH (F ROH=0.055) for Modicana, whereas Reggiana showed the lowest one (F ROH=0.035). The individuals of Italian Holstein and Reggiana showed high number of short ROH segments. Modicana and Cinisara showed similar results between them with the total length of ROH characterized by the presence of large segments due to a recent inbreeding. In this study, we reported the distribution of ROH islands across the genome of these cattle breeds to provide insights into the mechanisms underlying genomic differences among them.
Genomic regions with high frequency in runs of homozygosity
In our study, we did not perform LD pruning, but, owing to the minimum 4 Mb size of ROH segments, we tried to avoid small autozygous segments caused by LD. Indeed, a strong LD, typically extending up to ~200 kb, is common throughout the bovine genome (Mastrangelo et al., Reference Mastrangelo, Saura, Tolone, Salces-Ortiz, Di Gerlando, Bertolini, Fontanesi, Sardina, Serrano and Portolano2014), and short ROH are very prevalent. To exclude these short and very common ROH, the minimum length for ROH was set to >4 Mb.
The top 1% of SNPs with the highest number of occurrence was chosen as an indication of a possible ROH island in the genome. The same threshold was reported in studies on cattle (Szmatola et al., Reference Szmatola, Gurgul, Ropka-Molik, Jasielczuck, Zabek and Bugno-Poniewierska2016) and sheep (Purfield et al., Reference Purfield, McParland, Wall and Berry2017). Gaspa et al. (Reference Gaspa, Marras, Sorbolini, Ajmone-Marsan, Williams, Valentini, Dimauro and Macciotta2014) and Sölkner et al. (Reference Sölkner, Ferenčaković, Karimi, Perez O’Brien, Mészáros, Eaglen, Boison and Curik2014) used top regions with percentage of SNP in ROH >40% within breed, whereas Mészáros et al. (Reference Mészáros, Boison, Pérez O’Brien, Ferenčaković, Curik, Da Silva, Utsunomiya, Garcia and Sölkner2015) applied a threshold of 10%. Recently, a common ROH proportion higher than 7.5% was chosen as an indicator of potential autozygosity islands in cattle (Kukučková et al., Reference Kukučková, Moravčíková, Ferenčaković, Simčič, Mészáros, Sölkner, Trakovická, Kadlečík, Curik and Kasarda2017). Therefore, we have employed a stricter criteria compared with the last two works mentioned above.
The ROH peaks were distributed and shared among individuals, and it was clear that they were signs of common ROH islands within breeds. Some of these genomic regions overlapped with ROH islands found in other studies (Table 3). The ROH islands reported on BTA4 and BTA5 in Modicana overlapped with ROH islands reported in Pinzgau (Kukučková et al., Reference Kukučková, Moravčíková, Ferenčaković, Simčič, Mészáros, Sölkner, Trakovická, Kadlečík, Curik and Kasarda2017) and Simmental (Szamatola et al., Reference Szmatola, Gurgul, Ropka-Molik, Jasielczuck, Zabek and Bugno-Poniewierska2016). Several studies (Sölkner et al., Reference Sölkner, Ferenčaković, Karimi, Perez O’Brien, Mészáros, Eaglen, Boison and Curik2014; Mészáros et al., Reference Mészáros, Boison, Pérez O’Brien, Ferenčaković, Curik, Da Silva, Utsunomiya, Garcia and Sölkner2015; Szamatola et al., Reference Szmatola, Gurgul, Ropka-Molik, Jasielczuck, Zabek and Bugno-Poniewierska2016; Kukučková et al., Reference Kukučková, Moravčíková, Ferenčaković, Simčič, Mészáros, Sölkner, Trakovická, Kadlečík, Curik and Kasarda2017) showed ROH islands located on BTA6. These regions overlapped with the ones obtained in our study for Modicana (34.32 to 41.34 Mb) and Reggiana (38.69 to 39.35 Mb). Sölkner et al. (Reference Sölkner, Ferenčaković, Karimi, Perez O’Brien, Mészáros, Eaglen, Boison and Curik2014) studying Taurine and Indicine cattle breeds, identified a region in BTA16 (43.80 to 44.97 Mb) visible only in Taurine. This overlapped with a region obtained in our study in Cinisara (43.92 to 45.55 Mb). Similar results were also reported for Red Polish, Simmental, and Limousine cattle breeds (Szmatola et al., Reference Szmatola, Gurgul, Ropka-Molik, Jasielczuck, Zabek and Bugno-Poniewierska2016). Runs of homozygosity islands identified on BTA20 in Italian Holstein were also described by Szamatola et al. (Reference Szmatola, Gurgul, Ropka-Molik, Jasielczuck, Zabek and Bugno-Poniewierska2016) in Holstein. Moreover, the peak identified on BTA26 in Italian Holstein partially overlapped with those obtained for the same breed by Gaspa et al. (Reference Gaspa, Marras, Sorbolini, Ajmone-Marsan, Williams, Valentini, Dimauro and Macciotta2014) in which the stearoyl-CoA desaturase (SCD) locus is located. Finally, the ROH islands reported on BTA28 in Cinisara overlapped with an ROH island reported in Simmental (Szamatola et al., Reference Szmatola, Gurgul, Ropka-Molik, Jasielczuck, Zabek and Bugno-Poniewierska2016). These results suggested that some of the ROH islands are common among different cattle breeds, and harbor variants that are undergoing selection independently of production and selection characteristics (Szmatola et al., Reference Szmatola, Gurgul, Ropka-Molik, Jasielczuck, Zabek and Bugno-Poniewierska2016). The inconsistencies among the criteria defining ROH islands makes it difficult to compare studies because the lack of consensus allows different thresholds and thus different signals (Peripolli et al., Reference Peripolli, Munari, Silva, Lima, Irgang and Baldi2016). However, the overlapping ROH islands among studies provided good evidence that they are not artifacts but genuine genomic regions affected by inbreeding. It is important to highlight that ROH islands can also be partly explained by the reduced recombination rate. Indeed, despite ROH being more or less equally distributed throughout the chromosomes, ROH islands were mostly found in regions with low recombination rates (Purfield et al., Reference Purfield, McParland, Wall and Berry2017). To verify this distribution in our cattle breeds and to determine if recombination rate impacted ROH islands, the linkage information published by the USDA (Ma et al., Reference Ma, O’Connell, VanRaden, Shen, Padhi, Sun, Bickhart, Cole, Null, Liu, Da and Wiggans2015) were considered. Ma et al. (Reference Ma, O’Connell, VanRaden, Shen, Padhi, Sun, Bickhart, Cole, Null, Liu, Da and Wiggans2015) assessed the relationship between recombination rate and chromosomal locations because recombination rates are known to differ considerably across chromosomal locations, including telomeres and centromeres. All cattle autosomes are acrocentric with the centromere located at the beginning and the telomere at the end of the chromosome. These authors reported a very low recombination rate near the centromere and the beginning of each chromosome, and they showed that the middle of the chromosome had a decreased recombination rate, although the centromere is far. Moreover, Ma et al. (Reference Ma, O’Connell, VanRaden, Shen, Padhi, Sun, Bickhart, Cole, Null, Liu, Da and Wiggans2015) highlighted that this low recombination rate in the middle of chromosome was not universal across all bovine chromosomes, but more pronounced for some of them (i.e. BTA9, 10, 11, 13, 15, 16, 19 and 23). Following the smooth spline plotting of the recombination rate reported by Ma et al. (Reference Ma, O’Connell, VanRaden, Shen, Padhi, Sun, Bickhart, Cole, Null, Liu, Da and Wiggans2015) in cattle, we checked if the ROH islands shown in Table 1 overlapped with regions of the genome showing low recombination rates. Some ROH islands actually overlapped with regions having low recombination rates (Table 1) as reported in previous studies in sheep (Purfield et al., Reference Purfield, McParland, Wall and Berry2017). Moreover, a previous study on cattle (Purfield et al., Reference Purfield, Berry, McParland and Bradley2012) reported a correlation between extensive LD and high incidence of ROH. The majority of SNPs within ROH islands showed similar LD levels as those computed for the entire chromosome (Supplementary Material Table S2), with the exception of two ROH islands (on BTA16 in Cinisara and on BTA 6 in Reggiana). Therefore, their existence was not easy explained on the basis of just LD (Nothnagel et al., Reference Nothnagel, Lu, Kayser and Krawczak2009).
Comparison among overlapped runs of homozygosity (ROH) islands here detected and those reported in previous studies
Identification of candidate genes within runs of homozygosity
We found that some SNPs occurred in regions of poor gene content. Some of the identified ROH islands, such as on BTA10 in the Cinisara breed, contained only one annotated gene (WDR72) or uncharacterized genes (i.e. LOC107132862). This may reflect selection acting on uncharacterized regulatory regions or simply the fixation of non-coding DNA by genetic drift due to the absence of any selection (Qanbari et al., Reference Qanbari, Pausch, Jansen, Somel, Strom, Fries, Nielsen and Simianer2011). An enrichment of genes involved in several GO-terms was observed in the four cattle breeds. We have not discussed in detail all of the genomic regions associated with ROH islands. Instead, we focused on selected genes in highly GO-enriched terms that, on the basis of the literature, showed associations with several specific traits related to livestock. Therefore, the functions of candidate genes within ROH islands play important roles in cattle and other livestock species are summarized for each breed.
In Cinisara, the ROH islands were identified on BTA8, 10, 13, 16, 23 and 28. A total of 40 genes were identified as being related to catalytic activity (GO:0003824), with genes implicated in immune response and immune regulation (PIK3CD and SPSB1, respectively) (Ramey et al., Reference Ramey, Decker, McKay, Rolf, Schnabel and Taylor2013). A high number of genes (n=57) were identified as being related to cellular process (GO:0009987). Among these, some candidate genes mapped on BTA16, such as PEX14, which is related to dairy production, KIF1B, which is under strong selection in dairy Holstein cattle (Flori et al., Reference Flori, Fritz, Jaffrézic, Boussaha, Gut, Heath, Foulley and Gautier2009), and RERE, which is implicated in embryonic growth and reproductive development (Ramey et al., Reference Ramey, Decker, McKay, Rolf, Schnabel and Taylor2013). Moreover, 52 identified genes were also related to cell part (GO:0044464) in which we highlighted the ADK gene on BTA28, which is involved in a physiological state (Ramey et al., Reference Ramey, Decker, McKay, Rolf, Schnabel and Taylor2013). Other candidate genes within the ROH islands on BTA28 were NRG3 and PPYR1, which are related with bovine mammary gland development and milk production, respectively (Ogorevc et al., Reference Ogorevc, Kunej, Razpet and Dovc2009).
The ROH islands in the Modicana breed were identified on BTA1, 4, 5, 6, and 8. A total of 37, 82 and 69 genes were identified as being related to binding (GO:0005488), cellular process (GO:0009987) and cell part (GO:0044464), respectively. Several enriched GO-terms contained genes related with milk production, such as the LALBA gene, a major whey protein that showed a significant association with the milk protein profile (Huang et al., Reference Huang, Peñagaricano, Ahmad, Lucey, Weigel and Khatib2012). On BTA6, the most homozygous region (>45% of individuals having the ROH island) was found (6:37 019 972 to 39 069 719 bp) and it contained an intriguing element. A quantitative trait locus (QTL) on this chromosome affecting milk fat and protein concentrations has been reported (Zhang et al., Reference Zhang, Boichard, Hoeschele, Ernst, Eggen, Murkve, Pfister-Genskow, Witte, Grignola, Uimari, Thaller and Bishop1998). The QTL, containing six genes (ABCG2, PKD2, SPP1, MEPE, IBSP and LAP3), was identified within one ROH island in our study. In this chromosomal region, several genes associated with milk production traits are annotated, such as FAM13A1, a gene near a milk protein QTL related to the protein content (Cohen et al., Reference Cohen, Reichenstein, Everts-van der Wind, Heon-Lee, Shani, Lewin, Weller, Ron and Seroussi2014). The ABCG2 gene harbors a quantitative trait nucleotide for milk composition in cattle (Olsen et al., Reference Olsen, Meuwissen, Nilsen, Svendsen and Lien2008). LAP3 has been associated with milk production traits (Zheng et al., Reference Zheng, Ju, Wang, Li, Huang, Zhang, Zhong and Wang2011) and with calving ease in dairy cows (Olsen et al., Reference Olsen, Meuwissen, Nilsen, Svendsen and Lien2008). Other interesting candidate genes, based on their suggested molecular function, were found, such as CAV1 and CAV2 on BTA4, which are implicated in the immune system (Qanbari et al., Reference Qanbari, Pimentel, Tetens, Thaller, Lichtner, Sharifi and Simianer2014). On BTA5, the KRT gene family, which is associated with epithelial development, was highlighted, together with TFCP2, which contained a QTL associated with fertility (Moore et al., Reference Moore, Pryce, Hayes, Chamberlain, Kemper, Berry, McCabe, Cormican, Lonergan, Fair and Butler2016). Moreover, we observed a genomic region within the ROH island that contained olfactory receptor family genes. Olfactory receptors detect and identify a wide range of odors and chemosensory stimuli, a necessity for finding food, detecting mates and offspring, recognizing territories and avoiding danger. They are also reported to be duplicated within the bovine genome, suggesting that they may be under strong selection for newly evolving functions (Qanbari et al., Reference Qanbari, Pimentel, Tetens, Thaller, Lichtner, Sharifi and Simianer2014).
In the Reggiana breed, the ROH islands were identified on BTA1, 3, 6, 17, 26 and BTA29. A total of 30 genes were identified as being related to catalytic activity (GO:0003824) with candidate genes, such as DIRAS3, which is involved in reproductive traits (Cheng et al., Reference Cheng, Zhang, Deng, Jiang, Xiong, Li and Lei2007), and PTEN, which is involved in mammary gland function (Li et al., Reference Li, Xie, Wang, Bian, Li, Gao and Wang2015). A total of 44 genes were identified as being related to cellular process (GO:0009987) with genes, such as MINPP1, which is associated with milk fatty acid traits in dairy cattle (Li et al., Reference Li, Xie, Wang, Bian, Li, Gao and Wang2015). Moreover, 35 genes identified in ROH islands were related to cell part (GO:0044464), and several of these genes (such as TAOK3 and NCAPG) have been previously associated with milk production traits (Li et al., Reference Li, Wang, Moore, Schenkel and Stothard2010; Weikard et al., Reference Weikard, Widmann, Buitkamp, Emmerling and Kuehn2012). We also detected the SLC35D1 gene, which is associated with the immune system (Qanbari et al., Reference Qanbari, Pimentel, Tetens, Thaller, Lichtner, Sharifi and Simianer2014), and the KSR2 gene, which affects milk production traits (Pimentel et al., Reference Pimentel, Bauersachs, Tietze, Simianer, Tetens, Thaller, Reinhardt, Wolf and König2011).
In the Italian Holstein breed, the 11 ROH islands were found on BTA3, 5, 7, 10, 13, 20 and BTA26. The most representative supercluster identified by the GO-term enrichment contained several candidate genes involved in milk production traits and reproduction. Among the 86 genes related to catalytic activity (GO:0003824), we found RHOV, which is related to dairy production (Gutiérrez-Gil et al., Reference Gutiérrez-Gil, Arranz and Wiener2015), on BTA10 and SCD, which has a large influence on milk fat composition because it plays a major role in determining the monounsaturated fatty acids, primarily oleic acid and the CLA content of milk fat (Rincon et al., Reference Rincon, Islas-Trejo, Castillo, Bauman, German and Medrano2012), on BTA26. The greatest number of genes (138) was related to cellular process (GO:0009987). Of these, BMP4 is involved in the development and functioning of follicles (Qanbari et al., Reference Qanbari, Gianola, Hayes, Schenkel, Miller, Moore, Thaller and Simianer2010), and OXT and AVP play major roles in regulating estrous behavior in dairy cows (Kommadath et al., Reference Kommadath, Woelders, Beerda, Mulder, de Wit, Veerkamp, te Pas and Smits2011). Finally, 121 genes related to cell part (GO:0044464) were also found. The PELO gene on BTA20 is involved in dairy production (Gutiérrez-Gil et al., Reference Gutiérrez-Gil, Arranz and Wiener2015). Other important candidate genes within the ROH island were DIO1, which is related to milk synthesis and energy metabolism (Connor et al., Reference Connor, Sonstegard, Kahl, Bennett and Snelling2003) on BTA3 and SLC2A4RG, which is involved in lactation persistency (Nayeri et al., Reference Nayeri, Sargolzaei, Abo-Ismail, Miller, Schenkel, Moore and Stothard2016), on BTA10. C9, which is involved in immune response, was located within a QTL region for mastitis-related traits (Sahana et al., Reference Sahana, Guldbrandtsen, Thomsen and Lund2013).
As reported above, several enriched GO-terms were related to milk production, reproduction, immune response, and resistance/susceptibility to infections and diseases. This indicated that the analyzed individuals may have experienced selective pressure on their genomes for these specific traits. Some genomic regions may be fixed in individuals within a population as a result of artificial or natural selection for reasons such as adaptability or productivity. Cinisara and Modicana are two breeds that have excellent abilities to adapt to harsh environments, high resistance levels to infections and diseases, good maternal aptitudes, and high-quality milk production. Genes that are involved in these traits were detected in our study using the ROH approach and were consistent with the phenotypic characteristics of these two breeds. Recently, a study on local sheep breeds (Mastrangelo et al., Reference Mastrangelo, Tolone, Sardina, Sottile, Sutera, Di Gerlando and Portolano2017) revealed the presence of ROH islands in genomic regions that harbor candidate genes for selection in response to environmental stress and which underlie local adaptation. The presence of many immune system-related genes in the identified ROH islands could reflect selection (natural or artificial) for disease resistance. Reggiana and Italian Holstein are two breeds reared and selected for milk production, and in accordance with this phenotypic trait, our results emphasized the presence of dairy-related genes within the ROH islands. Currently, in dairy cattle, such as Holstein, the systemic decline of fertility is being observed, in agreement with the several genes implicated in affecting the reproductive traits highlighted in this work. Kim et al. (Reference Kim, Cole, Huson, Wiggans, Van Tassell, Crooker, Liu, Da and Sonstegard2013) found that several genomic regions within ROH were associated with economically important traits, including milk, fat and protein yields. Therefore, the annotated genes that mapped to these ROH islands were perceived as exposed to selection.
Conclusion
In this work, we examined the distributions of ROH islands across the genomes of four cattle breeds with similar production aptitudes but different selection histories. We confirmed that the ROH islands were clearly non-uniform among breeds and across chromosomes. In fact, different ROH islands were found across breeds, consistent with possible signatures of either artificial or natural selection. For most genes associated with ROH islands, a biological link to traits of economic importance, which are known to be under selection, can be hypothesized and are consistent with the phenotypic characteristics of these breeds. Because genomic regions that are subjected to selection tend to generate ROH islands, their distributions can indicate genomic regions that may have been subjected to selective pressure. Our results contributed to understanding how selection can shape the distribution of ROH islands and suggested that ROH islands can be used to identify genes potentially involved in economically important traits. Further research must be performed to compare selection signatures and ROH islands, and to incorporate the use of ROH island’ distributions across the genome to limit the number of false positives identified and to modify current procedures.
Acknowledgments
This research was financed by PON02_00451_3133441, CUP: B61C1200076005 funded by MIUR.
Declaration of interest
None.
Ethics statement
None.
Software and data repository resources
None.
Supplementary material
To view supplementary material for this article, please visit https://doi.org/10.1017/S1751731118000629
