The accumulation of harmful genes within the ROH hotspot regions of the Tibetan sheep genome does not lead to genetic load

Abstract

Background

Prolonged natural selection and artificial breeding have contributed to increased uniformity within the Tibetan sheep population, resulting in a reduction in genetic diversity and the establishment of selective signatures in the genome. This process has led to a loss of heterozygosity in specific genomic regions and the formation of Runs of Homozygosity (ROH). Current research on ROH predominantly focuses on inbreeding and the signals of selection; however, there is a paucity of investigation into the genetic load and selective pressures associated with ROH, both within these regions and beyond. On one hand, genes located situated ROH hotspot regions exhibit a degree of conservation in their genomic segments; on the other hand, these regions may also serve as critical loci for identifying signals of selection.

Results

High-throughput re-sequencing technology was utilized to investigate the ROH hotspot regions across 11 Tibetan sheep populations, resulting in the identification of ten conserved genes (ARHGEF16, Tom1l2, PRDM16, PEMT, SREBF1, Rasd1, Nt5m, MED9, FLCN, RAI1) that are associated with lipid metabolism, lactation, and development. These genes exhibited highly conserved within the ROH hotspot regions across all Tibetan sheep populations. Employing the integrated haplotype score (iHS) method, we screened for selective sites within frequently observed ROH hotspot regions to elucidate genomic differences among Tibetan sheep populations. A comprehensive analysis was conducted, involving Rnhom, dN/dS ratios, missense/synonymous ratios, and loss-of-function (LOF)/synonymous ratios, to investigate the accumulation of deleterious genes and the associated genetic load both within and outside ROH hotspot regions. The results revealed a higher accumulation of deleterious genes and a reduced genetic load within the ROH regions.

Conclusions

This study provides a comprehensive and precise genomic analysis and interpretation of Tibetan sheep, offering theoretical basis for genetic breeding and evolution in Tibetan sheep.

Introduction

Current research on the germplasm characteristics of Tibetan sheep has evolved from the analysis of phenotypic traits and physiological and biochemical indices to a genomic level, specifically emphasizing ROH. ROH are defined as long homozygous segments within the genome. The lengths of ROH segments in different individuals are correlated with the genetic background of their common ancestor; as the genetic relationship becomes more distant, the length of the ROH segments decreases [1]. ROH are influenced by various factors, including inbreeding, population bottlenecks, genomic characteristics, and artificial selection, all of which leaves distinct imprints on the genome [2, 3]. Throughout the sheep breeding process, selective breeding has targeted various traits, such as body shape, fur quality, meat quality, and milk production, leading to the emergence of numerous specialized sheep breeds. This selective pressure also affects the distribution pattern of ROH fragments within the sheep genome [4,5,6,7]. The substantial genetic information encapsulated within ROH offers a novel approach for investigating sheep genomic data. With the advancement in high-throughput re-sequencing technology, researchers have increasingly leveraged ROH to assess inbreeding, following the pioneering work of Ferencakovic et al. first employed this approach. This methodology is now frequently employed by scientists to evaluate the extent of inbreeding, infer population structure, and screen for functional genes associated with significant economic traits [8,9,10,11,12,13].

The genetic basis of inbreeding depression is rooted in the elevated levels of homozygosity observed in the genomes of inbred individuals, primarily arising from two mechanisms: recessive mutational influence and the loss of over-dominant contributions [14]. This phenomenon can diminish traits related to population adaptability [14, 15]. Specifically, Inbreeding negatively impacts not only adaptive traits such as survival [16], fecundity [17], and disease susceptibility [18], but also production traits, including milk yield and meat quality [19, 20]. Furthermore, a reduction in inbreeding can lead to alterations in serum ion concentrations and hormone levels, as well as changes in ATPase activity and various physiological and biochemical indicators [21, 22]. Inbreeding depression varies in terms of its timing, characteristics, and effects across different sexes, and its impact on species can differ significantly [14, 23, 24]. Notably, an intriguing study revealed that inbreeding,, when coupled with environmental interactions, increases the sensitivity of inbred individuals to oxidative stress and other environmental stressors, thereby heightening their risk of decline or mortality [25].

Genetic variation stemming from mutation, recombination, and gene flow enables organisms to consistently develop adaptive traits that enhance the fitness of their species [26, 27]. However, during this process, organisms may also fix deleterious mutations, leading to a reduction in species fitness. The decline in fitness in this context corresponds to the generation of genetic load [28]. The decline in fitness associated with inbreeding is partly due to genetic load, specifically inbreeding load, which is expressed exclusively in homozygous individuals; the hidden load arises from homozygosity, leading to increased exposure to deleterious alleles [29]. This inbreeding load is influenced by evolutionary factors such as purifying selection and genetic drift. Moreover, beneficial mutations can be introduced from outside the population through purifying selection and genetic rescue, which may alleviate the overall reduction in inbreeding load and enhance fitness [15]. Deleterious mutations include both missense mutations and loss-of-function (LOF) genetic variants, both of which are predicted to disrupt gene function and significantly impair an individual’s survival capacity. Consequently, accurate quantification of genomic missense mutations and LOF genetic variants is crucial for enhance our understanding of the impact of mutation load on species viability. In this study, we analyzed the genomic data of 11 Tibetan sheep populations to identify genes associated with economically significant traits. Our objective was to explore the differences in selection signals and inheritance patterns that arise from varying selection pressures during population formation. Additionally, we sought to understand the effects of genetic load on Tibetan sheep, providing a theoretical foundation for their genetic breeding.

Materials and methods

Sample collection, DNA extraction, and sequencing

A total of 220 blood samples were collected from 11 populations of Tibetan sheep across Gansu, Qinghai, and the Tibet Autonomous Region (Fig. 1A, Table S6). From each population, 20 adult Tibetan sheep were randomly selected, and 5 mL of blood was obtained for DNA extraction. The DNA concentration was measured, and the integrity of the DNA was evaluated using agarose gel electrophoresis. Genomic DNA was randomly fragmented into short segments using restriction enzymes, followed by end repair, the addition of dA tails, and ligation of sequencing adapters. The resulting DNA fragments were purified using AMPure XP beads, and segments ranging from 300 to 400 bp were selected for PCR amplification. Finally, the constructed library was purified and subjected to quality control before sequencing on the Hiseq X10 platform using a PE150 protocol, achieving a sequencing depth of approximately 5X.

Data filtering and genome alignment

The raw image data obtained from sequencing were transformed into sequence data, referred to as raw reads, through a process known as base calling, with the results stored in FASTQ file format. To ensure data quality, the raw reads underwent a quality control process, followed by rigorous filtering to obtain high-quality clean reads for subsequent analysis. Quality trimming was performed on short sequences, which primarily involved removing reads that contained adapters, reads with more than 10% N bases in single-end sequences, and low-quality reads in which over 50% of the bases had a quality score of Q ≤ 20.

After applying the aforementioned filtering criteria, the remaining high-quality read data were aligned to the Huoerba Tibetan sheep genome, which was accurately assembled by the Sheep Resources and Breeding Innovation Team at the Chinese Academy of Agricultural Sciences, using BWA software. The BAM files generated by BWA were then processed, with short sequence duplicates being removed using Picard software. To minimize mismatches in regions adjacent to InDels, local re-alignment was conducted using the Genome Analysis Toolkit (GATK) to obtain precise variant information, particularly for InDels. Additionally, GATK was employed to re-calibrate base quality, ensuring that the quality values of the read bases in the final output BAM files accurately reflected the true probabilities of mismatches with the reference genome, ultimately yielding high-quality and reliable variants. Genome coverage and sequencing depth were assessed using Bedtools (v2.27.1).

Mutation detection and annotation

The Unified Genotyper module of GATK (version 3.4–46) was employed to perform variant detection on the processed alignment files across multiple samples. The detected variants underwent filtering using the Variant Filtration tool with the parameters: -Window 4, -filter “QD < 4.0 || FS > 60.0 || MQ < 40.0”, and -G_filter “GQ < 20”. Here, QD denotes Variant Confidence/Quality by Depth; FS represents the Phred-scaled p-value derived from Fisher’s exact test for strand bias detection; MQ indicates Root Mean Square (RMS) Mapping Quality; and GQ refers to Genotype Quality. Finally, the functional consequences of the identified variants were annotated using ANNOVAR.

Analysis of genomic Ho (observed heterozygosity) and He (expected heterozygosity)

For the filtered SNP loci, the PLINK software was used to calculate the observed heterozygosity (H_o) and expected heterozygosity (H_e) of the population.

Detection of ROH

Prior to detecting ROH, indels were filtered from the markers using vcftools-0.1.14, and the original marker loci were filtered with PLINK. This study leveraged the genomic characteristics of different populations to conduct ROH detection in PLINK by employing a sliding window of specified SNP count across the chromosomes with a defined step size. Windows that met the specified criteria were then concatenated to complete the ROH detection.

The PLINK parameter settings included:

• 1. A minimum length of 500 Kb for detected ROHs;

• 2. A minimum SNP number of 20 for defining ROHs in each population;

• 3. A sliding window size of 10 Kb with an overlap ratio of 0.05;

• 4. A maximum allowed heterozygote count of 1 within the window;

• 5. A minimum SNP density of 1 SNP per 50 Kb within the ROH segment;

• 6. If the interval between two consecutive SNPs exceeds 50 Kb, they cannot be classified as part of the same ROH.

The results of the ROH detection were statistically analyzed across the different populations, focusing on metrics such as the number and length of ROHs, as well as F_ROH.

Formula of F_ROH:

$$\:{\text{F}}_{\text{R}\text{O}\text{H}}=\frac{\sum\:{\text{L}}_{\text{R}\text{O}\text{H}}}{{\text{L}}_{\text{A}\text{U}\text{T}\text{O}}}$$

$\:{\text{L}}_{\text{R}\text{O}\text{H}}$ is the length of ROH in the genome; $\:{\text{L}}_{\text{A}\text{U}\text{T}\text{O}}$ is the total length of autosomes covered by SNPs.

ROH hotspot region detection

The ROH segments detected in the 11 Tibetan sheep populations using PLINK were analyzed at the population level to calculate the ROH ratio (the proportion of SNPs located within ROHs) utilizing R software. Based on the calculated ROH ratios for the SNPs across the various populations, the top 1% of ROH ratios was designated as the threshold for high-frequency SNPs, with SNPs exceeding this threshold identified as ROH hotspot regions.ion of ROH hotspot regions.

Gene annotation of ROH hotspot regions

The ROH hotspot regions identified in the 11 Tibetan sheep populations were subjected to gene annotation utilizing the Ensembl database, with the gene annotation results documented in Table S1. Based on this annotated gene list, we conducted a targeted screening for genes associated with economic traits in Tibetan sheep, employing the Animal QTL Database, NCBI resources, and relevant published literature. Furthermore, the KOBAS online analysis tool was utilized to perform KEGG and GO analyses on the genes related to economic traits within the ROH hotspot regions.

Select signal detection

Using fastPHASE software, we performed haplotype phasing on the SNP data from the 11 populations of Tibetan sheep. The resulting output haplotypes were organized into the format required by the iHS software, utilizing the rehh package in R for the calculation of iHS scores. The computed iHS scores were standardized by subtracting the genome-wide average iHS from each score and then dividing by the standard deviation. Selection loci were identified by filtering for iHS values in the top and bottom 1% of the distribution, resulting in potential selectively favored loci that were subsequently annotated. Finally, an intersection was established between the identified genes and those located within the ROH hotspot regions, yielding the final set of selected genes.

Genomic Rnhom and dN/dS analysis

To evaluate the inbreeding depression status of a population through the calculation of the ratio of missense mutation inside and outside of regions of ROH. Rnhom was calculated by dividing the ratio of missense to synonymous counts of homozygous genotypes within ROH by the corresponding ratio outside ROH [30]. Additionally, the dN/dS ratios for ROH and non-ROH regions were computed using the following formula for dN/dS [31, 32]:

$$\:\text{d}\text{N}/\text{d}\text{S}=\text{d}\text{N}\left(\text{n}\right(\text{L}\text{O}\text{F})+\text{n}(\text{m}\text{i}\text{s}\text{s}\text{e}\text{n}\text{s}\text{e}\_\text{v}\text{a}\text{r}\text{i}\text{a}\text{n}\text{t}\left)\right)\div\text{d}\text{S}\left(\text{n}\right(\text{s}\text{y}\text{n}\text{o}\text{n}\text{y}\text{m}\text{o}\text{u}\text{s}\_\text{v}\text{a}\text{r}\text{i}\text{a}\text{n}\text{t}\left)\right)$$

Genetic load analysis

The Snpeff 4.3 software was utilized to annotates SNPs from 11 populations in relation to reference genomes, identifying three distinct types of mutation in regions of ROH and non-ROH: synonymous mutations, missense mutations, and LOF mutations. The genomic ratios of LOF to synonymous mutations and missense to synonymous mutations were calculated and subsequently mapped for both ROH and non-ROH regions [31, 33]. Furthermore, the genes harboring mutations were subjected to KEGG and GO enrichment analyses for additional insights.

Results

The total length of individual genomic ROH primarily falls within 1000 Mb, with the number of individual ROH totaling fewer than 1600, There is a linear relationship between the total length and the total quantity of ROH (Fig. 1C). The lengths of the ROH across the 11 Tibetan sheep populations increase with the number of SNPs. The correlation coefficient between the total number of ROH and the total length of ROH for individual Tibetan sheep genomes is high (Fig. 1C). The ROH length of the Tibetan sheep genome is mainly concentrated in 0–1 Mb, and the less ROH length is longer than 3 Mb (Table S7). In ROH regions, a higher number of SNPs indicates greater genetic variation and thus higher levels of genetic diversity. Notably, the ZSJ population exhibits a lower number of SNPs in its genomic ROH, while the GJ population displays a higher SNP count (Fig. 1B). The individual F_ROH values calculated from the genomic ROH can be used to assess the degree of inbreeding among individuals. It is evident that there is considerable heterogeneity among individuals in the Tibetan sheep populations. Based on F_ROH values, inbreeding levels can be categorized into three tiers: low (F_ROH < 0.1), medium (0.1 < F_ROH < 0.2), and high (F_ROH > 0.2). Individuals with F_ROH values greater than 0.2, indicating severe inbreeding, were identified in the GJ, OL, and GBB populations (Fig. 1D).

Fig. 1

Sampling and ROH Analysis Results of Tibetan Sheep Populations (A) Map of Tibetan Sheep Sampling Locations (B) Number of Homozygous SNPs in Tibetan Sheep Populations (C) Length and Quantity of ROH in Tibetan Sheep Populations (D) F_ROH Values of Individual Tibetan Sheep

ROH hotspot region genes in Tibetan sheep populations

This study investigates the distribution of ROH hotspot regions and their associated genes across 11 populations of Tibetan sheep, with findings summarized in Table S1. The number of identified ROH hotspot regions varies significantly among these populations. Specifically, the AW population exhibits the fewest ROH hotspots, totaling only 25, whereas the GBW population shows the highest number with 95 ROH hotspots. These ROH regions span multiple chromosomes and encompass a diverse array of genes.

The gene count within these ROH hotspot regions also differs among the populations, with the GBB having the highest number of associated genes, while the WT has the lowest (Fig. 2). Notably, the number of ROH hotspot regions within different chromosomes varies across the Tibetan sheep populations. Moreover, the gene counts within these ROH hotspots regions differ across the populations, with GBB having the highest number of associated genes, while WT has the lowest, as illustrated in Fig. 2. The distribution of ROH hotspot regions across different chromosomes also varies among the Tibetan sheep populations. The genes identified within these ROH hotspot regions are linked to various biological functions, including immune responses, lactation, reproduction, body weight regulation, and lipid metabolism. Notably, certain genes found in the ROH hotspots regions of Tibetan sheep populations adapted to higher altitudes—specifically the WT, ZSJ, GBW, and HB populations—are associated with mitochondrial function and adaptation mechanisms to high-altitude hypoxia.

Among these, several genes such as ARHGEF16, Tom1l2, PRDM16, PEMT, SREBF1, Rasd1, Nt5m, MED9, FLCN, and RAI1 are highly conserved within the ROH hotspot regions across Tibetan sheep populations. These genes are primarily involved in lipid metabolism, lactation, and developmental processes, underscoring their potential significance in the adaptive characteristics of these populations.

Fig. 2

Number of genes in ROH hotspot region of Tibetan sheep

Enrichment analysis results of ROH hotspot regions in 11 Tibetan sheep populations

We performed GO and KEGG enrichment analyses on the genes identified within the ROH hotspot regions of the 11 Tibetan sheep populations. The analyses were conducted with strict filtering criteria, including a minimum of two genes per enriched term and P-values of less than 0.05 for both KEGG and GO enrichment assessments. The results of the enrichment analyses reveal that the identified signaling pathways are predominantly associated with metabolic processes, molecular functions, disease pathways, and immune response mechanisms. Notably, the variation in the number of enriched gene entries across the different Tibetan sheep populations aligns with the varying gene counts detected within the respective ROH hotspot regions, as detailed in Table S2.

ROH related results of Tibetan sheep population

To assess the impact of inbreeding events on the accumulation of deleterious mutations, we calculated the Rnhom values, which represent the ratio of missense to synonymous mutations within ROH hotspot regions compared to the corresponding ratio outside of these ROH hotspot regions in the genomes of Tibetan sheep. The analysis revealed that the mean Rnhom values are significantly higher in the GJ, TS, and ZSJ populations, while the KC population exhibits the lowest mean value. Furthermore, the Rnhom values within these populations demonstrated considerable variability, particularly highlighted by the differences observed between the GJ and TS populations, as illustrated in Fig. 3A. These findings suggest that inbreeding events contribute to the accumulation of deleterious mutations in the genome. To quantify the selection pressure status of the Tibetan sheep population at the genomic level using non-synonymous mutation rate (dN) and synonymous mutation rate (dS), the dN/dS ratio can mitigate the effects of varying mutation rates across different genomic regions. In our study of 11 Tibetan sheep populations, the dN/dS ratios across all individuals were observed to be consistently less than one, indicating that these populations are under purifying selection. This suggests that the majority of non-synonymous mutations present in these populations are likely deleterious, with only a limited number being beneficial. Notably, the dN/dS ratios within ROH hotspot regions were found to be lower than those in non-ROH hotspot regions, implying that non-synonymous mutations occurring within ROH hotspot regions tend to be more detrimental than those outside these regions. Among the populations studied, both GJ and TS exhibited higher dN/dS ratios within ROH hotspot regions, indicating a relatively lower intensity of purifying selection compared to other populations, as depicted in Fig. 3B and C.

Fig. 3

Accumulation of Deleterious Mutations Load in the Tibetan Sheep Genome (A) Rnhom Values in the Tibetan Sheep Genome (B) Genomic dN/dS values of Tibetan sheep population (C) dN/dS values of ROH hotspot region in Tibetan sheep population

Inbreeding can increase the genetic load of a population’s genome to some extent, and genetic load can impair gene function, thereby reducing the population’s adaptability [34]. In the genomes of these Tibetan sheep, the genetic load within ROH hotspot regions was found to be lower than that in non-ROH hotspot regions, with missense/synonymous ratios being higher than LOF/synonymous ratios. This indicates that inbreeding in these Tibetan sheep populations does not necessarily lead to an increase in the genetic load of the genome, with a greater proportion of genetic load arising from missense mutations (Fig. 4A-B).

Fig. 4

Genomic genetic load of Tibetan sheep population. (A) The genome of Tibetan sheep is LOF/synonymous (B) The genome of Tibetan sheep is missense/synonymous

iHS analysis results for ROH hotspot regions in Tibetan sheep populations

To gain deeper insights into the selective pressures acting on loci within high-frequency ROH hotspot regions across various breeds of Tibetan sheep, we applied a haplotype-based iHS method to detect selection signals within these regions in 11 Tibetan sheep populations. The annotation results of the iHS analysis for the ROH hotspot regions are presented in Table S3, while Table S4 details the KEGG and GO analysis outcomes for the genes associated with the iHS signals within these regions.

We defined selection loci as the top 1% of SNP loci based on their iHS values from the Tibetan sheep populations, identifying candidate genes for the Tibetan sheep populations as those SNP loci shared by four or more groups, as summarized in Table 1 and illustrated in Fig. 5. The findings indicate that seven high-frequency ROH hotspot regions harbor SNPs under strong selective pressure, from which we identified nine candidate genes.

Table 1 Candidate genes of tibetan sheep population

Fig. 5

iHS analysis of loci within high-frequency ROH hotspot regions of the Tibetan sheep population

ROH gene mutations within the genome

The mutation load of genes associated with ROH in the Tibetan sheep genome was found to be enriched. The mutations across 11 Tibetan sheep populations were linked to various biological processes, encompassing metabolism, genetic information processing, cellular processes, organismal systems, and diseases pathways (Figure S1). A noteworthy observation from this study is that mutations predominantly affected organismal systems and diseases phenotypes, with a higher prevalence in genes implicated in cellular regulation. Specifically, the mutations in the COL20A1, MYO15B, METRN, and EHD1 genes were identified in over half of the Tibetan sheep populations (Table 2). These genes play crucial roles in the development, immune response mechanisms, and signal transduction in Tibetan sheep. COL20A1 is primarily found in bone and cartilage. It is an essential component for maintaining the structure and function of bone and cartilage [35]. METRN is predominantly expressed in the central nervous system (CNS) and is believed to play a significant role in immune-related regulation. EHD1 is a member of highly conserved gene family that encodes proteins containing the EPS15 homology (EH) domain. EHD1 is essential for the endocytosis of insulin-like growth factor 1 (IGF1) receptors and participates in various signal transduction pathways.

Table 2 LOF and missense mutations in the ROH hotspot region of Tibetan sheep

Discussion

The Tibetan sheep population generally exhibits inbreeding and low individual diversity

This study elucidated the current status of genomic inbreeding, genetic load, and selection pressure within Tibetan sheep populations, with the aim of providing valuable insights for future genetic breeding and conservation of biodiversity. The findings indicated a significant prevalence of inbreeding among Tibetan sheep, characterized by a high dispersion of F_ROH. Notably, populations such as GJ, OL, and GBB exhibited F_ROH values exceeding 0.2. A high inbreeding coefficient poses a risk of extinction and can lead to inbreeding depression [41, 42]. The inbreeding coefficient among certain individuals of Tibetan sheep is significantly higher than previously reported. Liu calculated the mean F_ROH for this population to be 0.0206, suggesting that the inbreeding status of Tibetan sheep has not improved and, in fact, has deteriorated over time [43]. Numerous studies have demonstrated that excessively high F_ROH values can have detrimental consequences, including an increased incidence of recessive diseases and a decline in population adaptability. For instance, in Holstein cows, elevated F_ROH values are correlated with reduced sperm concentration (SC) and sperm motility (SM), thereby significantly impairing reproductive performance [44]. Furthermore, in Arctic foxes, a 0.1 increase in F_ROH is linked to a 76% reduction in survival rates [45], while in wild Soay sheep, a similar increase in F_ROH corresponds to a 60% decrease in survival rates [46].

The reduction of heterozygosity caused by inbreeding will impact growth and reproduction

It is noteworthy that the genes identified in the ROH hotspot regions population of 11 Tibetan sheep are predominantly associated with body weight, lactation, reproduction, immunity, and metabolism. The majority of the enriched gene pathways are related to individual growth and development, as well as disease susceptibility and immune responses [47, 48]. A total of ten highly conserved shared genes were identified in these Tibetan sheep populations, namely Arhgef16, TOM1L2, PRDM16, PEMT, SREBF1, RASD1, NT5M, Med9, FLCN, and RAI1. These genes are implicated in a variety of biological processes, including cell morphology, transport, metabolism, transcription, milk production, and ontogeny [49, 50]. By collaborating with iHS to analyze the selection signals of the ROH hotspot regions locus, seven regions of Tibetan sheep were identified in which SNPs exhibited strong selection during the process of interrupted rearrangement. Nine genes (COMMD4, KCNT1, RNF217, YIPF1, TRIM67, NDC1, EBF4, CEP164, PTPRF) were detected in at least four populations. Notably, TRIM67, CEP164, KCNT1, and NDC1 are primarily associated with ontogeny. Research indicates that TRIM67 is involved in brain development, cytoskeleton regulation, and neuron morphogenesis, while CEP164 plays a role in embryonic development. KCNT1 is essential for regulating neuronal membrane excitability, and NDC1, a component of the nuclear pore complex, is significant in spermatogenesis [51,52,53,54,55,56]. PTPRF is integral to breast development, and its disruption has been associated with syndromic amastia [57]. The genes selected in this study are predominantly linked to ontogeny, genome stabilization, and the regulation of cellular metabolism. This information contributes significantly to our understanding the evolutionary trajectories of various Tibetan sheep populations.

Inbreeding maintains the conservation of adaptive genes to some extent

The ROH hotspot regions of Tibetan sheep has been specifically annotated with genes associated with immune system response and functional regulation (IRF2, NOD2, IL2, TLR9, CCL19, PIGR, ICAM2, CD19) [58,59,60,61,62,63,64,65], mitochondrial function (SLC25A11, MFN2, MTERF4, MIGA2, NDUFV1, TUFM) [66,67,68,69,70], and altitude hypoxia adaptation (COX4I2, EGLN1) [71, 72]. The conservation of these genes among Tibetan sheep can be attributed to inbreeding practices, coupled with the functional advantages conferred by these genes in high-altitude, low-oxygen environments. Individuals possessing these advantages genes are likely to experience enhanced survival and reproductive success, which can be understood as a consequence of natural selection. From an ecological and evolutionary standpoint, the preservation of these genes highlights the adaptive evolutionary processes of Tibetan sheep within high-altitude ecosystems. The selective pressure exerted on specific genes at elevated altitudes promotes their retention, thereby equipping populations to effectively address survival challenges. While inbreeding facilitates the maintenance of particular adaptive genes, it may concurrently result in diminished genetic diversity [14]. Low genetic diversity diminishes the adaptability of Tibetan sheep populations to environmental changes and may elevate the risk of recessive genetic disorder [73]. Therefore, breeding strategies for Tibetan sheep should be carefully evaluated to establish a balance between genetic health and overall fitness.

There was a higher accumulation of deleterious genes and a lower genetic load in the ROH hotspot region

This study revealed that the accumulation of deleterious mutations within regions of ROH was more pronounced compared to non-ROH hotspot regions of the genome. This finding suggests that inbreeding events contribute to the buildup of harmful mutations within the genome. In the Tibetan sheep genomes analyzed, the genetic load in ROH hotspot regions was found to be lower than that in non-ROH hotspot regions. Furthermore, the ratio of missense to synonymous mutations was higher than that of LOF to synonymous mutations. These results indicate that genomic inbreeding in Tibetan sheep does not necessarily lead to a rapid increase in overall genetic load, which appears to be primarily driven by missense mutations rather than LOF mutations. Relatedly, Kuang reported a similar genetic load phenomenon in Tibetan sheep while conducting research on endangered snub-nosed monkeys [74]. This perspective challenges the prevailing belief that increased inbreeding invariably results in negative genetic consequences, a notion that can be scientifically elucidated. The prompt removal of deleterious alleles during population bottlenecks in small populations can effectively mitigate mutation load [75, 76]. Van estimated the average deleterious nature of derived alleles across various mammalian species and found that those with historically small population sizes and low genetic diversity generally exhibited lower genetic loads compared to species with larger population sizes. Additionally, they observed that the process of genetic clearance evolves at a slower pace and that the accumulation of inbreeding mutation load occurs more gradually [74, 77]. These studies also support the notion that mutation load accumulates over a longer period than inbreeding [78].

Mutation load mainly affects organism immunity and organismal systems

Tibetan sheep predominantly inhabit high-altitude regions, where environmental stressors significantly enhance their immune adaptability. Research has demonstrated that variations in immune defense and stress levels exist among populations residing at differing altitudes, with those at higher elevations exhibiting a pronounced enhancement in innate immunity [79]. Animals with diminished immune function are unable to adapt to the harsh conditions of the Tibetan Plateau. Consequently, mutations or LOF in immune-related genes within the Tibetan sheep genome can directly affect their viability. The METRN mutation has been identified in more than half of the Tibetan sheep population. METRN is primarily expressed in the central nervous system and functions as a neuroprotective factor [80]. In addition, METRN is involved in the regulation of immune-related pathways, such as signaling through the B-cell receptor (BCR) and immunoregulation between lymphocytes [36]. Mutations in domestic animals residing at high altitudes can predispose them to altitude-related diseases. Recently, mutations associated with high-altitude pulmonary hypertension have been identified in beef cattle [81]. Furthermore, specific gene mutations have been detected in patients with β-thalassemia [82]. High altitude-related diseases can impact various systems, including the cardiovascular system, central nervous system, and reproductive systems [83,84,85]. These findings suggest that the decline in the Tibetan sheep population, attributed to inbreeding, primarily arises from the disruption of individual development and the functional integrity of immune-related genes, ultimately reducing their adaptability.