Genotyping data from 353 individuals across 18 Welsh breeds on the Illumina OvineSNP50 SNP array from Beynon et al. (2015) were used in this study. Illumina pair-ended read (150 bp) resequencing of 11 Welsh sheep samples from the same set, representing one sample per breed (Supplementary Table 1), was performed at the University of Aberystwyth to ∼13× raw coverage using Illumina HiSeq according to the manufacturer’s protocol. Three remaining Welsh sheep resequenced genomes for Welsh Hardy-Speckled Faced, Dolgellau Welsh Mountain, and Talybont Welsh Mountain were downloaded from the National Centre for Biotechnology Information Sequence Read Archive, PRJNA160933 (Heaton et al., 2014). A dataset of resequenced Russian sheep samples (n = 40) adapted to a contrasting environment was used as an outgroup in this study (Sweet-Jones et al., 2021).

Reads were mapped by using the Burrows-Wheeler Aligner BWA-MEM (BWA V.0.7.10 (Li, 2013) to the reference sheep genome Oar_V.3.1. Reads were sorted using the Samtools V.0.1.18 (Li et al., 2009) and duplicates were marked with the Picard V.2.18¹. Libraries were also merged using Picard. Base Quality Score Recalibration (BQSR) was performed, which account for the systematic errors in sequencing, using the Genome Analysis Toolkit (GATK V.3.8; McKenna et al., 2010). Samples then underwent variant calling with the GATK HaplotypeCaller function. Finally, all samples (n = 54) underwent joint calling to merge all reported variants into a single vcf file. Following this, hard filtering for quality scores assigned by BQSR was performed by GATK, using the filter expression “[QD < 2.0| | FS > 60.00| | MQ < 40.00| | MQRankSum < −12.5| | ReadPosRankSum < −8.0].” All variants from resequenced data were converted into a Plink V.1.90 (Purcell et al., 2007) format to be run through the DCMS pipeline (Yurchenko et al., 2019).

Welsh breed genotyping data were separated into three groups of related breeds based on the clustering analysis performed by Beynon et al. (2015), resulting in one group of upland breeds and two groups of lowland breeds. Of the two lowland groups, one consisted of five lowland breeds and the other consisted of only the Kerry Hill and Beulah breeds (KHB) (Table 1). Plink-formatted (Purcell et al., 2007) files for the upland and lowland breed groups had genotypes from the sex chromosomes removed and were filtered to remove the rare alleles [–maf 0.05], low called SNPs [–geno 0.01], or poorly genotyped individuals [–mind 0.05]. FastPhase V.1.4 (Scheet and Stephens, 2006) was used to estimate the number of haplotype clusters (k) for each group (Lowland k = 48, KHB k = 25, and Upland k = 53). HapFLK software (Bonhomme et al., 2010) was used to obtain selection statistics for each group. This test uses the hierarchal population structure to identify the haplotype-based selective sweeps, which focuses on the inherited combinations of alleles. It has the advantage of an increased statistical power, reliably detecting the hard and soft selective sweeps, and is a realistic simulation of selection through the haplotypes. HapFLK p-values were calculated using the Python script scaling_chi2_hapflk.py (Bonhomme et al., 2010; Fariello et al., 2013). Adjusted p-values, or q-values, were calculated through the R qqman q-value function (Turner, 2014). Selected intervals were determined by boundaries of q < 0.05 with SNPs within an interval of q < 0.01 considered to be under a strong selection.

Five established measures of selection and genetic diversity were both used on the genotyping (15 breeds) and whole genome resequencing (upland breeds n = 7; lowland breeds = 7) datasets as a DCMS (Table 1; Ma et al., 2015). Statistics used were: (i/ii) H2/H1 and H12, which can distinguish the hard and soft selective sweeps by measuring the intensity of selection (Garud et al., 2015); (iii) Tajima’s D comparing pairwise sequence differences and the number of segregating sites, detecting positive, negative, or balancing selection (Tajima, 1989); (iv) nucleotide diversity (Pi), average number of nucleotide differences between two sequences (Nei and Li, 1979); (v) F_ST fixation index, comparing single SNP frequencies across a population (Weir and Cockerham, 1984). By weighting the result of each statistic and generating a combined score, regions that overlap in the analysis outcomes gain a stronger evidence to be a region under selection.

Autosomal SNPs were filtered for –maf 0.0000001 –geno 0.1 –mind 0.1 by Plink. SNPs were then phased by chromosome using ShapeIt2 (Delaneau et al., 2011) with 400 states and an effective population size of 100. Phased chromosomes were split into appropriate groupings per chromosome using Plink and H2/H1, H12 were calculated using the H1_H12.py Python script (Garud et al., 2015). H2/H1 and H12 values were calculated in windows of 25 SNPs using a step size of one SNP following our previous study (Yurchenko et al., 2019).

Tajima’s D for mutation index was calculated over the same intervals of 25 SNPs, whose lengths were calculated from the output of the H2/H1 and H12 statistics. Using the vcftools V.0.1.13 (Danecek et al., 2011), Tajima’s D was calculated [–TajimaD 900000000] for each chromosome per group/breed file per window. Output files were concatenated per group.

F_ST was calculated with Plink comparing each group to all others. All negative values were converted to zero, and data were smoothed with the R runmed function in windows of 31 SNPs.

Plink-format file was split by chromosome per breed, and nucleotide diversity was calculated with the vcftools [–site-pi] option. Data were then smoothed using the R runmed function in windows of 31 SNPs.

Output files from individual statistics were sorted and joined by SNP id. Genome-wide p-values through ranking results of each statistic were calculated in the R MINOTAUR stat-to-p-value function specifying the one-tailed tests (H2/H1, H12, F_ST–right tailed, Pi, and Tajima’s D–left tailed). Covariance matrix was constructed based on sampling 300,000 randomly sampled SNPs using the R CovNAMcd function where the alpha = 0.75. DCMS statistics were calculated using the DCMS function and fitted to a normal distribution to examine normality implemented by the R MASS rlm function. These fitted DCMS values were converted to p-values using pnorm, and adjusted for a false discovery rate with the qvalue function. Q-values were parsed to determine selection with region boundaries set at a q < 0.2 and a threshold of q < 0.01.

For the regions defined by our pipelines as being under selection, genes were identified using a list of 26,958 genes of the Oar_V.3.1 genome downloaded from Ensembl BioMart v.98. Within each selected region, genes were then ranked based on their distances to the most significant SNP of that region with the closest gene being the top-ranking. The top 10 highest ranking genes from each region underwent literature review for their previous associations to adaptation to local environments or socioeconomic traits in animals. Genes with established links to these traits were identified as candidate genes in this study.

In the resequencing study, all SNPs were annotated with NGS-SNP (Grant et al., 2011) to identify missense mutations. Only genes with missense mutations in the regions under selection were considered for literature review. Additionally, missense SNPs with a strong support from the F_ST statistic (F_ST ≥ 0.3) were analyzed with PolyPhen2 to predict their effects on protein structure and function (0 = benign and 1 = deleterious; Adzhubei et al., 2013).

Functional enrichment analysis was performed using the Database for Annotation, Visualization, and Integrated Discovery (DAVID) v. 6.8 (Huang et al., 2009) using the same gene list downloaded from Ensembl Biomart v. 98. Enrichments were detected using the DAVID functional clustering tool, which verifies enrichments of similar GO terms across many different databases, to confer a stronger evidence of enrichment. Scores > 1.3-fold enrichment were considered.

CNV detection was performed to identify the regions in the genomes of Welsh sheep that had been duplicated or deleted with respect to the Texel reference genome. We identified CNVs in the resequenced samples with the cn.Mops R package (Klambauer et al., 2012) in a window length of 700 SNPs, using sequences that had undergone BQSR. This resulted in each individual being given a raw copy number (CN) per window. CN1–CN3 were considered to be normal and discounted from the results. Raw CNVs were merged into the CNV Regions (CNVRs) with the BedOps bedmap function using at least 50% reciprocal overlap in at least three individuals within the same group as criteria for inclusion. Duplicate CNVRs were removed, and the neighboring CNVRs were merged. This allowed us to be confident in our results by excluding the regions where CNVRs appear to overlap the signatures of selection, as these cast doubt over their reliability.

HapFLK clusterplots and haplotype trees were visualized from the prepared R scripts available (Bonhomme et al., 2010; Fariello et al., 2013). All Manhattan plots were rendered in R by the qqman package setting the suggestive line (q = 0.05) and genome-wide line (q = 0.01) to indicate selection. Haplostrips for visualization of haplotype sharing (Marnetto and Huerta-Sánchez, 2017) in the regions under selection detected by DCMS were run using phased data.

Using the HapFLK software, signatures of selection were found in each grouping using 44,711 SNPs (lowland group contained one region, KHB-two, and upland-31; Supplementary Tables 2–4). Cluster plots and Haplotype Trees for the most significantly selected regions in the KBH group are shown in Figure 1 and for other groups, in Supplementary Figure 1. Lengths of the regions found under selection ranged from 1.81 to 64.78 Mb with a median length of 7.33 Mb. Moreover, DCMS also found the regions under selection in all 15 breeds, with at least one region of a strong selection (q < 0.01) in each breed using 47,366 SNPs. Lengths of the regions under selection ranged from 1 to 3.9 Mb where the median length of regions was 0.16 Mb. The number of regions detected in each breed ranged from 8 in Kerry Hill to 89 in Clun Forest (Supplementary Tables 5–19). In these regions under selection, 1,089 unique genes were found, with 179 occurring in multiple breeds including: GBA3 (Hill Radnor, LWF, Lleyn, SWWM, TWM, and WMHF), ENSOARG00000021104 (BFWM, Clun Forest, Hill Radnor, and Lleyn), PCDH9 (BHC, Hill Radnor, and Llanwenog), and PPARGCA1 (LWF, SWWM, and WMHF). Manhattan plots for all the breeds investigated are shown in Figure 2 and Supplementary Figure 2.

There was a substantial overlap between the results for a region found on OAR6: 31.99–46.00 Mb (q = 2.0 × 10^–5) in the upland breeds in the HapFLK study to a narrower region, OAR6: 40.20–43.38 Mb, in LWF (q = 6.9 × 10^–15), WMHF (q = 4.3 × 10^–11), and SWWM (q = 8.73 × 10^–5). In all cases, the top ranked genes were GBA3 linked to liver metabolism (Dekker et al., 2011) and PPARGC1A with known roles in mitochondrial biogenesis, fat deposition, and milk fatty acid composition in cattle (Fernandez-Marcos and Auwerx, 2011; Armstrong et al., 2018; Yuan et al., 2019; Li et al., 2016). On the HapFLK cluster plot for this region, low diversity was seen in the Balwen, BWM, LWF, and SWWM breeds, further supporting the idea of selection for this region. When plotted on a haplotype tree, it can be seen that Balwen had the longest branch compared to the other upland breeds, whilst the strongest signature of selection is seen in SWWM (Supplementary Figure 1).

Several other candidate regions detected by HapFLK overlapped with the regions found in the Welsh breeds by DCMS. OAR3: 23.23–33.86 Mb (q = 5.0 × 10^–6) shared its top-ranked genes, TDRD15 and APOB, with a region under selection in the upland Balwen breed (q = 0.0003). These genes are associated with cholesterol mobilization in Large White pigs (Bovo et al., 2019). OAR7: 47.55–54.20 Mb (q = 0.0008) overlapped with another region OAR7: 51.30–52.76 Mb in the upland BWM (q = 5.0 × 10^–9), but these did not share the top-ranked genes. Strong overlapping candidate genes are presented in Table 2.

The most strongly selected region in the HapFLK study was seen on OAR15 32.15–72.10 Mb in the KHB lowland group, with the top-ranking gene being MUC15 (q = 4.1 × 10^–7), associated with a low fecal egg count in Spanish Churra sheep during gastrointestinal parasite infections (Periasamy et al., 2014; Benavides et al., 2015). This was supported by a low haplotype diversity seen in the Kerry Hill breed at this locus, but this signature was not seen in Beulah or in the DCMS results in either breeds.

In the other lowland group, the only region found under selection, OAR13: 47.16–54.37 Mb (q = 0.01), overlapped with another region found in the upland group, OAR13: 46.62–72.94 Mb (q = 2.0 × 10^–6), but the top candidate genes were different. In the lowland breeds, the top-ranked genes were RNF24 and PANK2, which have been found under selection in world sheep breeds in association to a loss of vision following domestication (Naval-Sanchez et al., 2018; Wang et al., 2019). Cluster plots for this region showed a decreased haplotype diversity amongst the selected region on OAR13 in the lowland breeds, especially Clun Forest and HSF. Furthermore, significant selection at this locus in HSF (q = 0.0006) was confirmed using the DCMS pipeline. In the case of the upland breeds, the top-ranked genes were DHX35 and PPP1R16B which are related to innate viral immunity (Rahman et al., 2017) and endothelial cell proliferation (Pszczola et al., 2018), respectively.

The genes found within the regions under selection from HapFLK, functional enrichments seen for the KHB group included the DENN domain and connexin gap junctions, enriched 1.6- and 1.4-fold, respectively. The most highly enriched terms in the upland breeds were bactericidal permeability protein, major intrinsic protein, and semaphorins, which were all enriched over threefold. The most enriched cluster seen in the genotyping DCMS analysis was the Type II keratin filaments from the Hill Radnor breed, which was enriched fivefold. This was followed by the Ribonuclease A in Clun forest, enriched 2.8-fold and leucine rich-repeats in Lleyn, enriched 2.7-fold (Supplementary Table 20).

Fifty-four (14 Welsh and 40 Russian sheep) resequenced genomes were aligned to the Oar_V.3.1 genome with a mean filtered coverage of 11.9× (Supplementary Table 1). A total of 41,643,098 SNPs were called, which were pruned to 38,276,494 SNPs after filtering. CNVRs covered 0.27% of the lowland and 0.24% of the upland genomes, overlapping 852 and 669 genes, respectively (Supplementary Tables 21,22). Three CNVRs from the lowland breeds and 52 CNVRs from the upland breeds overlapped the regions of selection. Some of these CNVRs had a high frequency in the population, including the regions under a strong selection on OAR24 in the lowland breeds, spanning the CLCN7 gene and on OAR17, whereas in the upland breeds, these spanned the IGLV4-69, ZNF280B, and PRAME genes.

After excluding the regions overlapping CNVRs, DCMS found 2,996 regions under selection in the Welsh breeds (lowland = 514, upland = 2,482; Supplementary Tables 23,24 and Figure 3A). These regions overlapped 104 and 430 genes in the lowland and upland breeds, respectively. The most significantly selected region in the upland group was OAR22: 15.3676–15.3679 Mb, which overlapped the NOC3L gene (q = 1.8 × 10^–8), and for the lowland group, it was a single synonymous SNP located in MYH11 (q = 0.0006).

In regions under selection, there were 12 missense mutations found in the lowland breeds (Supplementary Table 25) and 85 missense mutations found in the upland breeds (Supplementary Table 26). Of these, only missense mutations and their enclosing genes, which were top-ranking SNPs in their selected intervals are discussed below, leading to a total of 4 in the lowland breeds and 14 in the upland breeds. Clarification of the type of selection relied on a strong support from the H2/H1 and H12 statistics, which were considered to be haplotype-based selection, or F_ST, which was considered fixation of a variant in a population and was seen as a selection acting on individual SNPs.

Three missense mutations were found in the regions under selection (q < 0.01), with a strong support from the F_ST statistic (F_ST ≥ 0.3). In the lowland breeds, these included the reference alleles of OAR2:2,887,916 in the CDK5RAP2 gene (q = 0.006; F_ST = 0.3) and OAR4:47,087,846 in CDHR3 (q = 0.007; F_ST = 0.3). For the upland breeds, only one missense mutation was found, relating to the derived allele of OAR8:21,197,663 in the RWDD1 gene (q = 0.002; F_ST = 0.5). The Polyphen score of this selected mutation, L23P, was low and did not support the large change in the protein function (Figure 3B).

Missense SNPs in the regions under selection supported by the haplotype statistics (H2/H1, H12) were only found in the MC1R gene in the lowland breeds, and in the upland breeds: FGD3, HSPA5, ABCA13, PTPN1, ACP2, LOC101121718, NOC3L, VASH1, SIRPA, and UBA52 genes. Alternatively, some missense SNPs received strong support from Tajima’s D and Pi, including one in the GIGYF1 gene found in the selected regions in both the upland and lowland breeds, as well as the SLC22A7 and HSD3B7 genes in the upland breeds (Figure 3C).

In the upland breeds, three of the genes found with missense SNPs were also identified in the regions under selection defined in the upland breeds of the genotyping analysis. These included: VASH1 (q = 0.003) on OAR7, which also overlapped with BWM from the genotyping DCMS data; UBA52 (q = 0.01) on OAR13, which also overlapped with BHC from the genotyping data, and SIRPA (q = 0.007). No missense mutations found in the lowland breeds overlapped with the regions under selection detected from the genotyping dataset.

MC1R, found in a region under selection in the lowland breeds (Table 3), is known to cause an upregulation of tyrosinase in hair melanocytes, which led to an increase of black eumelanin pigment; however, in sheep, typically, white pheomelanin is produced due to selection for mutations in MC1R (Weatherhead and Logan, 1981; Yang et al., 2013). Another gene, CDK5RAP2, with a strong support from allele frequency statistics, is related to human neurodevelopment by recruiting tubulin subunits, which are important for cortical gyration (Issa et al., 2013). Mutations in this gene have also been linked to Hertwig’s anemia mutant mouse models displaying blood cytopenia, aneuploidy due to impaired cell-cycle spindle checkpoints, and increased neuronal cell death (Lizarraga et al., 2010). Finally, CDHR3, encoding a cell adhesion protein is linked to childhood asthma and rhinovirus-C susceptibility, was also located in a region under selection in Asian sheep by HapFLK (Fariello et al., 2014; Bochkov et al., 2015).

Within the upland breeds, many genes found in the selected regions had functions related to cell stress and metabolism (Table 3). The most highly supported by DCMS, RWDD1, encodes a transcription factor related to sulfide metabolism in metozoans (Kang et al., 2008; Li et al., 2018). HSPA5 also responds to cell stress mechanisms by promoting protein refolding in the endoplasmic reticulum in cancers or virally-infected cells (Booth et al., 2015). This gene has also been found under selection in Chinese Yellow-Feathered chickens with regard to meat quality (Huang et al., 2020) and muscling in world pig breeds (Li et al., 2011).

Of the three genes found in the selected regions in both the genotyping and resequencing datasets, two related to cell stress mechanisms. The most significant region contains VASH1, encoding the vasohibin 1 signaling molecule, with known roles of negatively regulating angiogenesis (Chen et al., 2020), as well as promoting expression of antioxidation enzymes (Miyashita et al., 2012). Age-related downregulation of VASH1 leads to a lower endothelial cell stress tolerance, posing as a risk factor for human vascular diseases in later life (Takeda et al., 2016). Secondly, UBA52 encodes a protein with an ubiquitinate activity, with its downregulation causing cell cycle arrest and reduced protein synthesis essential for pre-implantation embryogenesis success in mouse models (Kobayashi et al., 2016).

Several genes related to metabolism and growth were found under selection in the upland breeds. The most significant of these, HSD3B7, is part of the bile biosynthesis pathway (Cheng et al., 2003). PTPN1 is a risk-factor gene linked to diabetes and obesity (Olivier et al., 2004), which has a direct involvement in the insulin and leptin signaling pathways, and that mice lacking this gene were resistant to weight gain and intolerant to glucose (Elchebly et al., 1999). GIGYF1 also has roles in enhancing the insulin receptor pathway, but additionally has been linked to translational repression (Giovannone et al., 2003; Tollenaere et al., 2019). Similar effects have been seen with APC2, which is linked to muscle mass in mice (Kärst et al., 2011). The Rho-GEF-containing gene FGD3, expressed in the growth plate of long bones, was previously found under selection in French Trotter and Gidran horses, as well as in Jutland and Japanese black cattle in association to birth weight (Takasuga et al., 2015; Grilz-Seger et al., 2019; Stronen et al., 2019).

Two membrane transport proteins, ABCA13 and SLC22A7, were found to be in the regions under selection in the upland breeds. ABCA13 encodes a member of the ATP-binding cassette membrane transporter family, responsible for the active transport of biological substrates across cell membranes (Prades et al., 2002). Secondly, SLC22A7, a transmembrane solute carrier, with roles in cAMP and cGMP transport in mammalian tissues and, therefore, is important for intracellular signaling which may mobilize intracellular Ca²⁺, activate protein kinases, or activate transcription factors (Kobayashi et al., 2005; Yan et al., 2016). The final gene shared between the genotyping and resequencing data was SIRPA, which is expressed by macrophages and polarizes M1 phagocytic macrophages to M2 antiphagocytic macrophages, which is a key survival strategy for tumors (Barclay and van den Berg, 2014).

Nine functional category enrichments were found from genes within CNVRs in the lowland breeds and 14 enrichments in the upland breeds (Supplementary Tables 27,28). Some of these enrichments were shared between the lowland and upland breeds. These were semaphorins; ion transport, pleckstrin homology domain; SAND domains; and EGF-like domains. Exclusive enrichments in the lowland breeds’ CNVRs included cell surface receptors, neuromuscular process, and DNA binding whereas exclusive enrichments in the upland breeds included Src homology-3 domain Rho signal transduction, Notch signaling, and chondrocyte differentiation (1.4-fold). Genes within the regions under selection in the upland breeds showed significantly enriched clusters including interleukin-1 and ATP-binding (Supplementary Table 29). No functional category enrichments were found in genes in the regions under selection in the lowland breeds.