Comparison Among Approaches for Selection Signatures Detection
In this study, we investigated selection signatures in various Merino (and Merino-derived) sheep breeds using two different approaches. While the “FST-outlier” is considered a classical method for identification of regions putatively under differential selection in pairs of breeds (or group of breeds), the analysis of local ancestry, i.e., the genetic ancestry of an individual at a particular chromosomal location, in admixed populations to detect genomic regions where a significant excess of ancestry from a given parental breed exists (also known as “admixture mapping”) is so far a less popular approach. Among the “admixture mapping” approaches, LAMP has some interesting features that prompted us to opt for this method. Unlike algorithms that are based on reference haplotype frequencies for each of the parental populations, for which larger sample sizes are required to capture haplotypic diversity, LAMP relies on reference allele frequencies (Shriner, 2013) and is consequently less affected by a reduced sample size. Also, it operates on sliding windows of contiguous SNPs, using a “majority vote” for each locus, over all windows that overlap with the SNP, in order to decide the most likely ancestral population at the marker. This simple approach has been shown to provide fast and robust estimates (Sankararaman et al., 2008). Despite LAMP has been developed for estimation of the locus-specific ancestry in recently admixed populations, it has been shown to be robust to inaccuracies in the parameter “number of generations since the admixture.” A critical issue, limiting the widespread use of LAMP, is represented by the choice of the external reference samples to be used as proxies for the true ancestral populations, as the latter are generally not available for sampling (Shriner, 2013). In this study, a set of four hypotheses, each including a test breed and two proxies for the parental populations, were formulated based on historical knowledge on the origin of breeds and the inferred proportions of global admixture. These have to be interpreted with caution given the possible influence of complex patterns of historical admixture known among Merino and Merino-derived sheep populations (Ciani et al., 2015). The four hypotheses were hence tested using the algorithm implemented in LAMP. On the other side, for the “FST-outlier” approach, we were able to define, a set of nine pair-wise comparisons by contrasting (i) Merino populations of Iberian origin (Spanish Merino and Australian Merino, respectively) with non-Merino populations of Iberian origin (Churra and Ojalada, respectively), (ii) Merino populations of Italian origin (Gentile di Puglia and Sopravissana, respectively) with non-Merino populations of Italian origin (Appenninica and Bergamasca, respectively), and (iii) a Merino population of Asian origin (Chinese Merino) with a non-Merino population of Asian origin (Tibetan). The rationale behind the above pairing is that differentially selected loci may be easier to detect when contrasting more homogeneous breeds, such as Merino versus non-Merino breeds from the closest geographical area (Manunza et al., 2016).
Consistently with expectations, the two adopted approaches produced only partly overlapping signals. Indeed, the two methods rely on different algorithms and different assumptions, which also imposed a different organization of the dataset used with the two approaches (four single-breed tests vs. nine pair-wise comparisons, for the “local ancestry” and the FST-outlier” approaches, respectively), thus hampering direct head-to-head comparisons. Notwithstanding, the majority of the regions detected using the “local ancestry” method showed overlapping signals with the “FST-outlier” results. Although we interpreted the above as evidence supporting the robustness of the obtained results, it must be taken into consideration that regions identified by “local ancestry” were generally large, and significant SNPs detected via the “FST-outlier” method may likely occur in there by chance. Indeed, the number of loci identified as putatively under selection pressure using the “FST-outlier” method largely exceeds the number of putatively selected regions identified using the “local ancestry” approach.
In this study, we also investigated genomic regions of high homozygosity (ROH islands), as these have been shown to be abundant in regions under positive selection (Metzger et al., 2015; Mastrangelo et al., 2017; Purfield et al., 2017; Talenti et al., 2017a; Talenti et al., 2017b; Mastrangelo et al., 2018). While we observed little overlapping between ROH islands and regions identified via “local ancestry,” some overlapping was observed between ROH islands and SNPs detected as significant using the “FST-outlier” approach. ROH islands may be the consequence of the genetic hitchhiking phenomenon at loci physically linked to the variant site under direct positive selection pressure. The “local ancestry” approach looks for regions with an excess of ancestry from one of the two parental populations, and not necessarily these regions have to display high homozygosity, although this feature is likely to be observed in case of a strong selective sweep. Similarly, in the “FST-outlier” approach, homozygous genotypes (for different alleles in the two breeds) at loci physically linked to the variant site under direct positive selection pressure may display high frequencies if a strong differential selection existed among the two considered breeds. Also, the argumentation reported above may apply here: the larger number of loci identified as putatively under selection pressure using the “FST-outlier” method may be more likely to occur by chance within large ROH islands compared to the fewer genomic regions identified via “local ancestry.” Moreover, ROH analysis might detect selection related to any trait, while contrasting Merino and non-Merino is more likely to detect signals related to this specific trait. Finally, the existence of ROH islands may be due to non-genetic factors such as demography.
Best Candidate Regions and Putatively Selected Genes
As the aim of this study was to identify loci most likely associated with the Merino phenotype, we arbitrarily identified (i) the best candidate regions detected via the “local ancestry” approach as those being shared by all the four breeds and (ii) the best candidate SNPs detected via the “FST-outlier” approach as those observed in at least 70% of the pair-wise comparisons (six out of nine). Based on the above, two large regions on OAR17 and OAR18 were retained for (i), and three, on OAR3, OAR10, and OAR13, for (ii). The robustness of the adopted procedure was also suggested by the occurrence, in all of the five regions, of QTLs/associations known to be related to wool traits in the ovine species. Moreover, at OAR17, combining analysis of “local ancestry” and inspection of the sheep QTL database allowed to significantly shorten the candidate interval. On the contrary, known QTLs for wool traits described on OAR18 and OAR3 are mapped within extremely large chromosome intervals. These were identified by Allain et al. (2006) and Ponz et al. (2001) who performed whole-genome scans using microsatellite markers on experimental flocks obtained crossing Lacaune with Sarda, and Berrichon du Cher (a Merino-derived breed) with Romanov (a non-Merino breed), respectively. In what follows, the gene content of the best candidate regions is presented by chromosome order.
OAR3
The region detected on OAR3 contains five genes, LOC105609945 (long noncoding RNA), MSRB3 (methionine sulfoxide reductase B3), LOC105609947 (long noncoding RNA), LOC105609948 (a pseudo-gene), and LEMD3 (LEM domain containing 3). Interestingly, the sub-region containing the genes MSRB3, LOC105609947, and LEMD3 was found to harbor a selection signature putatively for tail fat deposition in previous studies contrasting thin- vs. fat-tail sheep breeds, from China (Yuan et al., 2017), and from North Africa and the Chios island (Mastrangelo et al., 2019), for adaptation when contrasting the Red Maasai sheep with the Ethiopian Menz (Fariello et al., 2014), and for ear morphology in Chinese sheep breeds (Wei et al., 2015) and in French Suffolk sheep (Rochus et al., 2018). The latter suggest to consider the genes encoded by the signal on OAR3, notably MSRB3 and LEMD3, as candidates for ear size based on literature showing the possible role of the two genes in ear position in dogs (Vaysse et al., 2011) and ear size in pigs (Wilkinson et al., 2013). A more detailed discussion of each single gene in the OAR3 selection signature is provided below.
LOC105609945—No evidence for involvement of LOC105609945 in any peculiar Merino feature was found.
MSRB3—The methionine sulfoxide reductase B3 (MSRB3, alias DFNB74) catalyzes the reduction of free and protein-bound methionine sulfoxide to methionine. This antioxidant repair enzyme has been described in human epidermal keratinocytes and melanocytes, as well as in hair follicles (Taungjaruwinai et al., 2009). It has been shown to be expressed also in inner and outer hair cells of mouse inner ear (Ahmed et al., 2011). Diseases associated with MSRB3 include deafness (https://www.genecards.org). Down-regulation of MSRB3 has been shown to impair the normal auditory system development through hair cell apoptosis in zebrafish (Shen et al., 2015). The gene has been found in previous selection signatures studies in sheep (Kijas et al., 2012; Fariello et al., 2014; Wei et al., 2015; Manunza et al., 2016; Yuan et al., 2017; Rochus et al., 2018; Mastrangelo et al., 2019). More interestingly, for this study, it has been found within a selection signature observed contrasting fine-wool Merino and coarse-wool Churra sheep breeds (Gutiérrez-Gil et al., 2017). Another line of evidence for the involvement of MSRB3 in hair/wool physiology comes from the observation that actin’s polymerization properties and actin cytoskeletal-mediated events, such as correct bristle development, which are altered by specific oxidation of its conserved methionine (Met)-44 residue on the pointed-end of actin subunits, are rescued by a methionine sulfoxide enzyme reductase (SelR/MsrB) in Drosophila (Hung et al., 2013). In this species, actin plays a role not only in bristle but also in wing hair development (Guild et al., 2005). In mammals, actin has been shown to be one of the major components of both the water-soluble and -insoluble fraction from hair and hair follicles (Vermorken et al., 1981; Lee et al., 2006). Actin bundles in the hair follicle would act as stress fibers and serve as a tensile scaffold for the growth and integrity of the hair follicle (Furumura and Ishikawa, 1996). In Tibetan sheep, microRNAs differentially expressed in wool follicles during anagen, catagen, and telogen phases, thus potentially regulating wool follicle development, targeted, among others, genes in the pathways that regulate the actin cytoskeleton (Liu et al., 2013). MSRB3 also contained the (intronic) SNP that, in this study, displayed the highest number of “FST-outlier” pair-wise comparisons showing signal of selection observed in the OAR3 region.
LOC105609947—No evidence for involvement of LOC105609945 in any peculiar Merino feature was found.
LOC105609948—It’s a pseudo for the ubiquitin-conjugating enzyme E2 D3 gene (UBE2D3), which is part of the bone morphogenic protein (BMP) signaling pathways (gene ontology database accession ID: GO:0030509). BMP ligands (BMP2 and BMP4) when expressed in dermal macro-environment during telogen (resting phase of hair cycle) have been shown to strongly suppress ability of resting hair follicles to be reactivated and grow again (International Patent no. WO2010059861A1 available at https://patents.google.com/patent/WO2010059861A1).
LEMD3—As previously mentioned, the LEM domain containing three gene (alias MAN1) has been found in various selection signatures studies in sheep (Fariello et al., 2014; Wei et al., 2015; Manunza et al., 2016; Yuan et al., 2017; Rochus et al., 2018; Mastrangelo et al., 2019), including the study by Gutierrez-Gil et al. (2017) where fine-wool Merino and coarse-wool Churra sheep breeds were contrasted. Moreover, it has been found associated with the abnormal hair quantity phenotype from the HPO Gene-Disease Associations dataset (Köhler et al., 2014).
OAR10
The region detected on OAR10 includes four genes: LOC106991357 (long noncoding RNA), LOC101110773 (elongation factor 1-alpha 1-like), RXFP2 (relaxin/insulin-like family peptide receptor 2), and LOC106991379 (a pseudo-gene). Despite this region was detected as putatively selected in studies investigating tail fat deposition (Moioli et al., 2015; Yuan et al., 2017; Mastrangelo et al., 2019) and adaptation (Yang et al., 2016; Seroussi et al., 2017), RXFP2 is the most studied gene and is well known for being involved in horn presence/absence and morphology in sheep (Johnston et al., 2011; Kijas et al., 2012; Fariello et al., 2014; Pan et al., 2018). A genome-wide association study for wool traits in Chinese Merino sheep detected, on this chromosome, a significant SNP for “fiber diameter” at position 30 Mb, together with several SNPs significant for “crimp” at 26–27 Mb (Wang et al., 2014). In what follows, a more detailed discussion of the four genes annotated in the OAR10 region is provided.
LOC106991357—No evidence for involvement of this locus in any peculiar Merino feature was found.
LOC101110773—It codes for an elongation factor 1-alpha 1-like. The elongation factor 1-alpha 1 (EF1A1) is a GTP-binding protein which has a primary function as an essential house-keeping gene by delivering aminoacyl-tRNAs to the ribosome during the elongation step of protein translation. EF1A1, together with genes associated to the Usher syndrome, a congenital disease characterized by perturbation of normal organization and growth of hair bundles within the inner ear, is a downstream target of GBX2, which induces EF1A1 activation upon binding to the EF1A1 core promoter. GBX2 has been shown to be expressed in the otic placode, which develops into the inner ear. Loss-of-function and mis-expression studies have shown that GBX2 is essential for development of the inner ear sensory organs. However, neither direct evidence for involvement of EF1A1 in hair bundles organization and growth nor in any peculiar Merino phenotype has been found so far. Another elongation factor type (EF1Bγ) has been proposed to bind to keratin (Kim et al., 2006). The presence of large amounts of EF1Bγ in keratin bundle rich hair fibers would support its biological role in the intermediate filament organization (Sasikumar et al., 2012).
RXFP2—This gene is involved, among others, in the biological process “activation of adenylate cyclase activity” (https://www.uniprot.org/uniprot/Q8WXD0). Adenylate cyclase is responsible for the synthesis of 3’,5’-cyclic adenosine monophosphate (cAMP). Agents that increase cAMP levels have been shown to be potent inhibitors of human and mouse hair follicle growth (Harmon and Nevins, 1997). However, we cannot exclude that, in our tested scenarios, different alleles at this gene may have been differentially selected in the considered breeds as a consequence of selection toward different horn phenotypes. Rochus et al. (2018) highlighted that a number of single nucleotide polymorphisms exist in French sheep in the region extending 100 Kb upstream of RXFP2, with haplotypes in polled sheep being distinct from those observed in horned sheep. From these findings, they suggested that multiple ancient mutations, rather than a single mutation, are likely affecting horn phenotypes. Pan et al. (2018) reported strong association between three SNPs within the RXFP2 gene and horn sizes in a Tibetan population characterized by the presence of animals with heterogeneous horn types.
LOC106991379—No evidence for involvement of this locus in any peculiar Merino feature was found.
It is worth mentioning that, only 0.076 Mb upstream to LOC106991357 on OAR10, the gene FRY is mapped (interval 28,959,450–29,212,913 bp). Looking at our results separately for each tested scenario, we observed that this interval was overlapping with the region detected by the “local ancestry” method in Chinese Merino, as well as with the regions detected by the “FST-outlier” method in the Chinese Merino vs. Tibetan, Australian Merino vs. Churra, and Gentile di Puglia vs. Appenninica comparisons. Moreover, the FRY interval was only slightly upstream to the regions detected by the “FST-outlier” method in the Australian Merino vs. Ojalada (0.2 Mb), Spanish Merino vs. Ojalada (0.18 Mb), Gentile di Puglia vs. Bergamasca (0.18 Mb), Sopravissana vs. Appenninica (0.2 Mb), and Sopravissana vs. Bergamasca (0.18 Mb). In sheep, FRY has been suggested as a key candidate gene for the piebald phenotype in Merino (Garcia-Gamez et al., 2011) and has been suggested to be associated with the black spot phenotype in Valley-type Tibetan sheep (Wei et al., 2015), and with differences in coat color pigmentation distribution between the Awassi and Afec-Assaf sheep (Seroussi et al., 2017). On the contrary, Zhang et al. (2013) detected FRY when contrasting Rambouillet and Suffolk sheep and suggested it to be a candidate gene affecting wool quality. Indeed, FRY encodes a protein furry homolog that, in Drosophila, has been found in growing hairs (He et al., 2005), and whose disruption has been shown to provoke abnormally branched bristles and strong multiple-hair phenotype, with clusters of epidermal hairs and branched hairs (Cong et al., 2001). Fang et al. (2010), following the transgenic FRY protein in vivo, found it to be highly mobile and to accumulate at the distal tip of growing bristles and suggest that it could function in directing/mediating the intracellular transport needed for polarized growth.
OAR13
The region detected on OAR13 contains a single gene (RALY). It encodes a heterogeneous nuclear ribonucleoprotein that binds poly-U-rich elements within several RNAs and regulates the expression of specific transcripts (Cornella et al., 2017; Rossi et al., 2017). In early 90s, the gene was shown to be involved in the pleiotropic lethal yellow phenotype of the mouse due to a deletion of the genes RALY and EIF2S2 (eukaryotic translation initiation factor 2 subunit 2), upstream the ASIP (agouti signaling protein) gene, responsible for the ectopic over-expression of the agouti signaling protein under the control of the RALY promoter (Michaud et al., 1993; Duhl et al., 1994). The gene has been repeatedly detected, often together with the neighbor ASIP gene, in association studies concerning skin pigmentation and skin neoplasms, (https://www.ncbi.nlm.nih.gov/gap/phegeni?tab=1&gene=22913; Jacobs et al., 2015), as well as in several other type of cancers, where it is considered to represent a metastatic marker (Roberts et al., 2019). In 2013, a mutation in this gene has been associated with the saddle tan and black-and-tan phenotypes in Basset Hounds and Pembroke Welsh Corgis (Dreger et al., 2013). Similarly, it has been associated with coat color phenotypes in Chinese and Iranian goats (Guo et al., 2018; Nazari- Ghadikolaei et al., 2018). Worth mentioning that a SNP at OAR13 (position 62.9 Mb) was associated by Bolormaa et al. (2017) with wool fiber diameter in Merino sheep. Although the SNP is not far from the RALY gene, it mapped within the EIF2S2 gene, which has been shown to be involved in protection against chemotherapy-induced alopecia (Nasr et al., 2013).
OAR17
Out of the two pairs of genes flanking the two QTLs on OAR17, one (LOC101120019) is a pseudo-gene related to a 60S ribosomal protein (L10A), one (TMEM132B) codes for a trans-membrane protein, one (LOC101115905) codes for refilin-A (alias FAM101A), and one (LOC106991703) is responsible for the production of a long noncoding RNA. Their possible involvement in the Merino phenotype is discussed here based on evidences from the literature.
LOC101120019—A mutation in a 60S ribosomal protein (L21) has been shown to be involved in hereditary hypotrichosis simplex (HHS), a form of nonsyndromic inherited hair loss disorders (Zhou et al., 2011). 60S ribosomal proteins (L6 and L24) have been shown to be expressed in human anagen hair samples (Carlson et al., 2018). Interestingly, the 60S ribosomal protein L10A has been shown to be expressed in root hairs of Medicago truncatula (Covitz et al., 1998). As is common for genes encoding ribosomal proteins, multiple processed pseudo-genes of the 60S ribosomal protein L10A are dispersed through the genome (https://www.ncbi.nlm.nih.gov/gene/4736). In particular, LOC101120019 on OAR17, being a pseudo-gene, is more likely to play, if any, a regulatory function on the hair physiology.
TMEM132B—The trans-membrane protein 132B is required for normal inner ear hair cell function (https://www.genecards.org/cgi-bin/carddisp.pl?gene=TMEM132E). TMEM132A, but not TMEM132B, TMEM132C, or TMEM132D, was found to be expressed in wool follicle bulb of Tibetan sheep during phase transformation from the middle anagen, to catagen and late telogen/early anagen (Liu et al., 2015). TMEM132E was found to be highly expressed in murine inner hair cells, and a variant in TMEM132E was identified as the most likely cause of autosomal-recessive nonsyndromic hearing loss. Knockdown of the TMEM132E ortholog in zebrafish affected the mechanotransduction of hair cells. (Li et al., 2015).
FAM101A—The gene product is involved in the regulation of the perinuclear actin network and nuclear shape through interaction with filamins. It plays an essential role in the formation of cartilaginous skeletal elements (UniProtKB:Q5SVD0). In addition, it has been shown to be differentially expressed in hair follicle stem cells residing in the bulge of mouse hair follicles versus the epithelial basal cells outside the bulge (Chang, 2014). FAM101A mRNA was detected via next-generation sequencing in wool follicle bulb samples of Tibetan sheep from middle anagen, catagen, and late telogen/early anagen phases (Liu et al., 2015). In a genome-wide association study performed using 50 K SNPs in a Baluchi sheep population, one of the significant SNP markers associated with greasy fleece weight was located within FAM101A (Ebrahimi et al., 2017).
LOC106991703—No evidence for involvement of LOC106991703 in any peculiar Merino feature was found.
OAR18
The region identified on OAR18 included a large number of genes (70) that obviously hampered a detailed analysis of the available literature for each single gene. We hence opted for checking which of the above genes could be retrieved by querying the human GeneCards database using keywords representative of Merino phenotypes. While none of the genes was retrieved when using “wool” or “horn” keywords, 16 genes were retrieved when using the keyword “hair” and, among them, three displayed particularly high GeneCards relevance scores (NFKBIA, SEC23A, and PAX9).
NFKBIA (NFKB inhibitor alpha)—It has been shown to modulate WNT, SHH, and LHX2IS signaling at early stages of hair follicle development in mice. In particular, in the epidermis of mice lacking the transcription factor nuclear factor-kappa B activity, primary hair follicle pre-placode formation is initiated without progression to proper placodes (Tomann et al., 2016). The gene has been also detected as putatively under selection in a Chinese Merino sheep population (Liu et al., 2017).
SEC23A—It is one of the major components and markers of COPII vesicles from endoplasmic reticulum. It has been found associated with the sparse hair phenotype in humans (https://mseqdr.org/hpo_browser.php?8070) Moreover, it may contain causative mutations for an autosomal recessive disease known as cranio–lenticulo–sutural dysplasia, alias Boyadjiev–Jabs syndrome, in which patients have abnormal hair, among other cranio-facial abnormalities. Also, it has been shown to co-localize with the three proteins, transmembrane (Cdh23), scaffold (harmonin), and actin-based motor (Myo7a), whose defect is responsible of various types of the Usher syndrome, a multi-genic congenital disease characterized by perturbation of normal organization and growth of hair bundles within the inner ear (Blanco-Sánchez et al., 2014)
PAX9—It is a member of the paired box (PAX) family of transcription factors. Heterozygous mutations in PAX9 have been associated in humans with non-syndromic tooth agenesis, non-syndromic, and familial oligodontia, with peg-shaped laterals and microdontia incisors. Often, these symptoms are associated with hair defects (Roberts and Chetty, 2018) as the same genes responsible for tooth development are involved in the growth and development of the other tissues derived from the ectoderm, including hair.
In general, biological evidence for the involvement in plausible Merino phenotypes was observed for the vast majority of coding genes in putatively selected regions detected either via “local” ancestry” or the FST-outlier” approach. The above result highlights the power of the multi-cohort approach adopted here. While we cannot exclude that false positive signals may have been retained in this study, still this represents so far the most complete genome-wide study of selection signatures for the Merino phenotype. The selection signatures reported here provide a comprehensive insight into the genetic basis underlining the Merino phenotype in sheep, which appeared here to be mainly represented by wool (and horn) traits. Targeted studies at both physiological and molecular levels will be needed to better understand the biological complexity behind these commercially relevant traits.