Loci Associated with Negative Heterosis for Viability and Meat Productivity in Interspecific Sheep Hybrids

Simple Summary One of the main aims of selection is enhancing the economically important traits of animals and plants. Crossbreeding is one of the common ways to improve the economically important traits. Usually, different species of one genus are used for crossbreeding. The use of different closely related species is known as hybridization. It is expected that hybrids will have advantages over their parents in terms of economic productivity due to the effect known as heterosis. However, the phenomenon of negative heterosis exists and can induce negative impacts on different important traits, such as viability, reproductivity, productivity, etc. The exact biological mechanisms of negative heterosis remain unknown. In this study, we aimed to determine the genetic factors associated with negative heterosis in interspecific hybrids between domestic sheep (Ovis aries) and argali (Ovis ammon). We discovered one novel locus associated with viability and two novel loci associated with meat productivity. For loci associated with meat productivity, we demonstrated the effect of negative heterosis. These results may help to understand one of the possible genetic mechanisms of negative heterosis. Abstract Negative heterosis can occur on different economically important traits, but the exact biological mechanisms of this phenomenon are still unknown. The present study focuses on determining the genetic factors associated with negative heterosis in interspecific hybrids between domestic sheep (Ovis aries) and argali (Ovis ammon). One locus (rs417431015) associated with viability and two loci (rs413302370, rs402808951) associated with meat productivity were identified. One gene (ARAP2) was prioritized for viability and three for meat productivity (PDE2A, ARAP1, and PCDH15). The loci associated with meat productivity were demonstrated to fit the overdominant inheritance model and could potentially be involved int negative heterosis mechanisms.


Introduction
Producing new breeds is an essential process for increasing the phenotypic indices of economically important traits, and hybridization is one of the ways to do so. Usually, different species of one genus are used for hybridization, but in some cases different closely related species are used. For example, in sheep genetics there are a few examples of hybridization of domestic sheep (Ovis aries) with their wild ancestors, such as the Asian mouflon (Ovis orientalis) [1], Pamir argali or Marco Polo sheep (Ovis ammon polii) [2], European mouflon (Ovis orientalis musimon) [3], and bighorn sheep (Ovis canadensis) [4]. It is expected that hybrids will have advantages over their parents in terms of economic productivity due to the effect known as heterosis, e.g., for the intraspecific hybrids of domestic sheep, the proportion of individually retained heterosis was the most important genetic factor associated with increased lamb survival [5]. Hybridization with wild ancestors is usually performed for increasing the body size, survivability, disease resistance, etc. For instance, hybrids of domestic sheep with argali showed increased lean meat production [1], while hybrids with bighorn sheep demonstrated resistance to Mannheimia haemolytica in the absence of vaccination [4].
Negative heterosis, on the other hand, can have negative effects on different important traits, such as viability, reproductivity, productivity, etc. Unfortunately, despite its economic importance, this phenomenon remains understudied.
In quantitative genetics there are several models of heterosis. The genes causing this phenomenon may affect traits in three different ways: dominance (attributes heterosis to the canceling of deleterious or inferior recessive alleles contributed by one parent, by beneficial or superior dominant alleles contributed by the other parent in the heterozygous genotypes at different loci) [6][7][8], overdominance (attributes heterosis to the superior fitness of heterozygous genotypes over homozygous genotypes at a single locus) [8,9], and epistasis (contribution of epistatic interactions between non-allelic genes) [10,11]. In the field of plant genetics, different studies favor one mechanism over another [12][13][14][15].
As for domestic animals, the majority of studies on negative heterosis and its possible genetic and cytogenetic mechanisms have been conducted on chickens [16,17]. These studies have found that non-additive genes, their related oxidative phosphorylation, and the difficulties in homology matching between the DNA sequences of genetically divergent breeds might lead to negative heterosis. As for other domestic animals, there are only several publications, but they do not provide sufficient information about the possible mechanisms of this phenomenon [18,19]. In any case, to claim that a given associated locus is involved in heterosis, it is not enough to simply prove that it is associated with the trait-the locus should also have a non-additive inheritance model with negative overdominance. Such a model is a classic example of negative heterosis, when heterozygotes show inferior performance compared to homozygotes.
Previously, our group focused on obtaining a population of hybrids between domestic sheep (Ovis aries) and argali (Ovis ammon) [20]. In this study, a strong negative heterosis was observed for interspecific hybrids relative to their diagonal body length and their chest width and depth. In addition, many hybrid lambs died in the first 3 months after birth, reflecting negative heterosis with respect to viability. Additionally, considering the fact that the genetic and molecular mechanisms of negative heterosis in sheep are unknown, the focus of the present study was on determining the genetic factors associated with negative heterosis in interspecific hybrids between domestic sheep (Ovis aries) and argali (Ovis ammon).

Animals
A total of 89 hybrids of F1 Romanovskaya × argali backcross were analyzed. Handling and breeding of the Russian sheep population followed the protocols approved by the Ethics Committee on Animal Care and Use of the Institute of Cytology and Genetics (approval # 45/2 of 10 January 2019).

Genotype Data
The animals were genotyped using the Ovine Infinium ® HD SNP BeadChip containing 603,350 (HD) SNPs (Illumina Inc., San Diego, CA, USA). The alleles were coded using the Illumina A/B scheme. Quality control of the hybrid genotypes was performed using information about four additional groups of animals (argali-10 animals; Katahdin-14; Romanovskaya-18; F1 hybrids of Romanovskaya and Katahdin-14). Earlier, some SNPs had been noted to have a very low quality of genotypes in each of these groups. These SNPs were removed from the analyzed hybrid genotypes to prevent methodological errors. Then, a standard quality control procedure was carried out, according to the following SNP and sample inclusion criteria: SNP and sample call rate > 0.98; sample heterozygosity within ±3 standard deviation ranges; and minor allele frequency (MAF) > 0.01. As a result, 87 animals and 446,790 SNPs passed the quality control.
To calculate a kinship matrix, a subset of the pruned SNPs that were in approximate linkage equilibrium with one another was generated using such options as maf 0.05-indeppairwise 500 kb 10.3 from PLINK 1.9 [21]. In total, 46,471 pruned SNPs were generated.
The kinship matrix was calculated using the IBs function from the R package Gen-ABEL [22].

Phenotype Data
Eleven carcass traits (see Supplementary Table S1) were measured four times-after birth, at 6, 42, 90, and 180 days, identifying 44 phenotypes in total. Based on these traits, four integral and holistic indices (see Supplementary Table S2) were calculated that related to meat productivity in sheep and cows [23,24]. These four indices and body mass were only used in relation to the 42-day timepoint as the latest one to reasonably affect the animals' survival rate. The missing phenotypes were imputed using the R package PHENIX [25]. For the imputation, all 11 original traits for 6 and 42 days after birth were applied.
If an animal had no phenotypic information starting from a given time, it was considered to be deceased by that timepoint. The numbers of living animals for each timepoint are presented in Table 1.

Genome-Wide Association Study (GWAS)
For both viability and meat productivity traits, an association test with two degrees of freedom was carried out, which did not restrict the locus inheritance model. The probable inheritance model was assessed by either association boxplots or by the lowest-association pvalue for one-degree association tests for additive, recessive, dominant, and overdominant models. The test was applicable only for the SNPs with all three genotypic classes available. If an SNP had only two genotypic classes, a one-degree-of-freedom test without inheritance model inference was applied.

Viability GWAS
We used the Cox regression model from the R package SURVIVAL. The likelihood ratio test (LRT) was applied for estimating association p-values based on the following procedure. The null model (H0) was set as follows: Survivability~sex + number of fetus + Genetic principal component 1 + Genetic principal component 2 (1) where survivability is the status of the animal (i.e., alive or dead) and when it died (see Section 2.3 for more details), and the 4 covariates (sex, number of fetuses, and genetic principal components 1 and 2) are based on the variance-standardized relationship matrix. For the alternative model (H1), SNP genotype was added as a factor. The association statistic was calculated as follows: where L( ) is the corresponding model's likelihood. The resulting statistics were distributed as chi-squared with one degree of freedom in the event that the SNP model had only two genotypic classes, or with two degrees of freedom in the event that it had all three genotypic classes. The p-value significance threshold was determined as follows: 0.05/Nindep snps = 1.08 × 10 −6 (3) where Nindep snps = 46,271 SNPs remaining after pruning. For a locus that passed the threshold, Cox regression using the mlreg function from the R package GenABEL was applied. This analysis was conducted for the four inheritance models (additive, dominant, recessive, and overdominant).

GWAS of the Four Indices and Mass
For this study, a mixed linear regression model was applied. The GRAMMAR-Gamma method was used to correct the phenotypes in relation to their relatedness structure [26]. To perform this correction, an identity-by-state (IBS) matrix was calculated as described above. Sex, number of fetuses, the last time the phenotype information was available, the two first genetic principal components, and SNPs associated with viability were considered as covariates. The LRT was applied to the resulting residuals as described above. The null model was calculated via linear regression as follows:

Phenotype~1
(4) and the alternative model was defined as follows: phenotype~SNP (5) where the SNP was determined as a factor. The p-value significance threshold value was set as follows: 0.05/(4 × Nindep snps) = 2.7 × 10 −7 (6) where Nindep snps = 46,271 SNPs remaining after pruning, and 4 is the number of independent traits used in the analysis.
The associated loci were defined as regions within ±500 kb around the SNPs most significantly associated with the trait (lead SNPs).
A locus was defined as new if it was not present on the list of known SNPs and genes from our previously published database of the QTLs and genes associated with meat productivity in sheep [27]. This database includes information from 24 papers published between 2013 and 2020.
Inflation factors for five traits are presented in Supplementary Table S3.

Functional Annotation
The Ensembl Variant Effect Predictor (VEP) was applied to analyze the potential effects of the SNPs and indels that were in high linkage disequilibrium (LD; r2 > 0.7) with the lead SNPs. The analysis was performed using the software available online [28].

Estimation of Heritability
For estimation of additive heritability, we used GREML analysis from GCTA software [29]. Sex, number of fetuses, the last time the phenotype information was available, the two first genetic principal components, and SNPs associated with viability were used as covariates.

Replication Analysis
For replication analysis, different types of hybrids of the Romanovskaya breed and Katahdin breed with argali and mouflon (35 animals in total) were used [30]. We performed the same protocol for GWAS as described above. The threshold for replication was set to 0.05/(4 × 2) = 0.00625 (7) where 4 is the number of independent traits and 2 is the number of discovered SNPs from GWAS of the four indices and mass.

Viability GWAS
The resulting inflation factor for the GWAS was 1.071. One locus significantly associated with viability ( Table 2) was found. The study's Manhattan plot is presented in Figure 1, and the regional association plot is shown in Supplementary Figure S1. Since viability is not a quantitative trait, the association boxplots were uninformative. For that reason, we could only speculate about the locus's inheritance model based on the association results for the four models investigated (additive, dominant, recessive, and overdominant). The lowest p-value (Supplementary Table S7) was observed for the additive model, with a p-value of 4.6 × 10 −3 , and the G-allele estimated effect was 0.051.

Functional Annotation
Detailed information on the functional annotation is presented in Supplementary Tables S5-S6. Two intron variants (rs417431015 and rs402808951) were discovered for the viability and compactness indices, respectively. The SNP of rs413302370 was located in the intergenic region.

GWAS of Four Indices and Mass
The GWAS's results are presented in Table 1. Two loci associated with the compactness index were detected; joint Manhattan plots for the compactness index and the boxplots for the two loci can be seen in Figure 2. Manhattan plots for four other traits are presented in Supplementary Figures S2-S5. Regional association plots for the two loci can be found in Supplementary Figures S6 and S7. As the boxplots show, both variants fit a strong overdominant inheritance model. The estimated additive heritability for the traits varied from 0.46 to 0.85, with the length index having the highest heritability (see Supplementary  Table S4).

Discussion
In the present study, three genetic variants associated with productivity traits in sheep hybrids were found, and four genes were prioritized in these loci. Locus rs417431015, associated with viability, had quite a significant effect on the given sample size. Using the variance approximation for the quantitative traits, written as Z 2 /N, we es-

Functional Annotation
Detailed information on the functional annotation is presented in Supplementary  Tables S5 and S6. Two intron variants (rs417431015 and rs402808951) were discovered for the viability and compactness indices, respectively. The SNP of rs413302370 was located in the intergenic region.
For the locus tagged by rs417431015, associated with viability, the nearest gene ARAP2 was prioritized. This gene is involved in a pathway that controls focal adhesion through other signal proteins [31]. It is also involved in the integrin α5β1 traffic pathway [32]. ARAP2 influences sphingolipid metabolism through glucosylceramide synthase [33]. Furthermore, it is associated with carcass traits in Santa Ines sheep [34], loin strength in Chinese Holstein cows [35], and carcass traits in Korean cattle [36].
As for the two loci associated with meat productivity (rs413302370 and rs402808951), three genes (PDE2A and ARAP1 for locus rs413302370 and PCDH15 for locus rs402808951) were prioritized. PDE2A encodes the phosphodiesterase 2A protein, and its mutations are associated with movement disorders in humans, such as chorea [37] and paroxysmal dyskinesia [38]. ARAP1 is involved in the regulation of membrane trafficking and actin cytoskeleton reorganization through activation of Arf and Rho GTPases [39]. In addition, ARAP1 is known as a common regulatory variant for type 2 diabetes in humans [40]. Mutations in the PCDH15 gene cause Usher syndrome type 1F in humans [41].

Discussion
In the present study, three genetic variants associated with productivity traits in sheep hybrids were found, and four genes were prioritized in these loci. Locus rs417431015, associated with viability, had quite a significant effect on the given sample size. Using the variance approximation for the quantitative traits, written as Z 2 /N, we estimated the viability variance for this locus to be about 29%. Since the most probable inheritance model for this locus was an additive one, it could not be involved in negative heterosis. For this locus, the ARAP2 gene involved in focal adhesion was prioritized. In addition, this gene is associated with meat productivity and carcass traits in cattle and sheep. Thus, additional investigation is required to clarify the mechanism of mutation in the ARAP2 gene, which may lead to death in interspecific sheep hybrids.
For the two loci associated with the compactness traits, a good fit to the strong overdominant inheritance model was observed. Their variances were relatively high as well, being 30.4% and 35.2% for rs413302370 and rs40280895, respectively. Interestingly, we found associations only for the compactness index, while three out four other traits demonstrated higher estimated additive heritability (length index, mass, and cumulative index). This may be related to the observed non-additive effects of SNPs associated with the compactness index and the phenomenon of negative heterosis. The overdominant model indicated the classic pattern of a genotype associated with a trait affected by heterosis, i.e., the heterozygotes showed substantially inferior performance compared to the homozygotes ( Figure 2B,C). The effective allele frequencies were significantly different between the parental breeds, which also supports the hypothesis that these loci were involved in negative heterosis. Moreover, these results confirm our previous observations of negative heterosis manifestations in relation to the diagonal body length and the chest depth and width, due to the associations demonstrated for the compactness index to comprise the diagonal body length and chest girth.
Unfortunately, there are only a few examples of studies of interspecific sheep hybrids in the literature. To the best of our knowledge, there have been no studies describing the phenomenon of negative heterosis in interspecific hybrids. Our results showing negative heterosis contradict the findings of a previous study that showed increased lean meat production for interspecific sheep hybrids with Asian mouflon [1], which can be explained by the fact that different domestic sheep breeds were used for crossbreeding in the experiments. Additionally, in our study, a different set of carcass traits was used. However, in another study where domestic sheep were crossed with European mouflon (Ovis orientalis musimon), the hybrids had less weight than domestic sheep [3]. Our results overlap with those of a study of interspecific sheep hybrids, where domestic sheep were crossed with Pamir argali or Marco Polo sheep (Ovis ammon polii) [2]. In that study, the authors prioritized PCDH10 as a candidate gene associated with body weight traits. In our study, we prioritized the PCDH15 gene for the compactness index, which belongs to the same protein family as PCDH10. Further investigations of these proteins are needed to understand the possible mechanisms of influence on meat productivity traits.
It should be noted that there have been a few sheep studies where heterosis was shown and measured. The effects of direct heterosis (i.e., due to a lamb being crossbred) were small for an Australian population including more than 1 million animals [42]. In another study, a large direct and maternal effect of heterosis (due to the lambs' dams being crossbred) was estimated for a U.S. sheep population [43]. In addition, the effect of positive heterosis was shown in a study crossing indigenous fat-tailed and commercial sheep breeds [44]. These studies show a possible positive heterosis in sheep and, in combination with our results, may help to make interspecific crosses between domestic sheep and their wild relatives more successful in terms of increasing their economically important traits.
It is noteworthy that for both viability and meat productivity the genes from one protein family (ARAP1 and ARAP2) were prioritized. Both of these genes are involved in membrane transport: ARAP1 regulates the cell-specific trafficking of a receptor protein involved in apoptosis, while ARAP2 encodes a protein that is associated with focal adhesions and functions downstream of RhoA to regulate focal adhesion dynamics. Presumably, there is a gene network including these genes, but additional investigation into interspecific sheep hybrids is required to clarify this issue.
F2 backcrosses are preferable to observe all three genotypic classes for each SNP; it is clear that the studied population did not have the most optimal design for heterosis investigation, since it consisted of F1 backcrosses and their mother breed (Romanovskaya). For that reason, we do not expect that our findings could be replicated using other populations, considering that the possible heterosis mechanisms could be unique for each hybrid population. We performed replication analysis using other types of hybrids (35 animals in total). No SNPs reached the significance threshold of replication.
Better understanding of the genetic causes of negative heterosis in different domestic animals may help to improve breeding results. Knowledge of "proper" heterosis-associated allele combinations may help us to develop efficient marker-assisted systems that can be applied for breeding purposes.

Conclusions
In summary, the findings of the present study suggest two possible genetic markers (with three prioritized genes) that could potentially be involved in negative heterosis mechanisms in interspecific hybrids between domestic sheep and argali. In addition, we discovered one novel locus that was associated with viability. Further investigation is needed to clarify the possible molecular role of the discovered genes in negative heterosis. The results of the present study can help researchers and breeders in futures studies to create new interspecific sheep hybrid populations that are more successful in terms of economically important traits, as well as to develop efficient marker-assisted systems for breeding purposes.
Supplementary Materials: The following supporting information can be downloaded at: https://www. mdpi.com/article/10.3390/ani13010184/s1, Figure S1: The regional association plot for rs417431015; Figure S2: The Manhattan plot for shortlegged index; Figure S3: The Manhattan plot for length index; Figure S4: The Manhattan plot for mass; Figure S5: The Manhattan plot for cumulative index; Figure S6: The regional association plot for rs413302370; Figure S7: The regional association plot for rs402808951;