Introduction

Evolutionary mechanisms of adaptation to environmental conditions are developed and genetically determined (1-3). These mechanisms are fundamental in extreme or near-extreme climatic conditions, such as hot climates typical of arid zones (4,5). Genetic selection is a crucial method to enhance productive traits in local sheep through crossbreeding and selection programs (6,7). The progressive DNA research methods made it possible to study local genetic resources, identify the features and patterns of economically valuable traits of animals at the molecular genetic level, and thereby move to a new stage of breeding work, considering the preservation and improvement of the breeding traits of the livestock. The methods have provided great potential for investigating genetic diversity and structure as well as unraveling the shared ancestral history of livestock populations (8-11). Sheep breeding in the Russian Federation is an important economic sector as a source of wool, meat, and dairy products (12). At present, the largest sheep populations in the country are located in the southern regions of Dagestan, Kalmykia, the Stavropol Territory, the Astrakhan Region, Karachay-Cherkessia, the Rostov Region, the Republic of Tyva and the Volgograd Region. The top ten sheep farms are closed by the Republic of Bashkortostan and the Republic of Altai. Body weight is one of the most essential and economically advantageous traits characterized by complex inheritance; therefore, searching for genetic mechanisms that affect its formation is of increased scientific interest (13).

Therefore, more attention should be given to the factors that affect productive efficiency, including genetic variations and the rearing and management of herds. In this case, the genetic structure in local herds of three sheep breeds, popular in the South of Russia, and its association with body weight were first investigated (GWAS) using a 50 K SNP chip. Considering the relevance, novelty, and significance of such studies, we aimed to assess the genetic diversity and genetic structure of three local sheep breeds in the South of Russia, including the Kalmyk fat tailed, Edilbay, and Volgograd breeds, as well as perform an association analysis intended to search for genes and genome regions associated with the body weight index.

Materials and methods

Ethical approve

The ethical approve to conduct this scientific work was given by Department of Animal Husbandry and Breeding work of the Committee of Agriculture of the Volgograd region and the Republic of Kalmykia (EA NIIMMP # 1-2022-01-10).

Animals and sampling

The Edilbay breed of meat-fat productivity was created on the territory of semi-desert and steppe pastures of the Urals and Volga at the end of the 19th century (14). There are still some unclarified issues regarding the breed's origin. Historical facts indicate that the Edilbay sheep were obtained from crossing local Kazakh fat-tailed sheep with rams of Kalmyk breed. Sheep highly adapted to the continental climate were selected for breeding. The advantages of this sheep breed are a strong constitution, a well-developed physique, a fat tail, and high precocity (12,15). In 100-day fattening, the average daily gain is 195 g, and the maximum is 253 g. Sheep have higher wool productivity than other fat-tailed sheep, i.e., the average wool shear is 3.5 kg in rams and 2.3-2.6 kg in ewes. Wool is characterized by heterogeneity and contains 52-56% fur fiber, 16-19% transitional hair, and 24-28% beard hair. The fineness of fur fiber reaches 18 microns, transitional hair – 33 microns, and beard hair – 60 microns. The fecundity of Edilbay sheep is low and amounts to 110-120%, with milk production being good and the livability of the young stock being high. Moreover, due to its precocity and good adaptive properties, the Edilbay breed is widespread throughout Russia. It is used to improve other breeds in terms of fattening and meat indicators and create new breed types of sheep (Figure 1).

Figure 1: Edilbay sheep breed.

In Kazakhstan, specialists from the South-Western Research Institute of Animal Husbandry and Crop Breeding have developed a meat-fat Ordabasin sheep breed by crossing ewes of the local Kazakh fat-tailed coarse-wooled breed with rams of the Edilbay and Hissar breeds (16). The Volgograd sheep breed originated in 1932 when crossbreed offspring were obtained by the method of the complex reproductive crossing of course-wool fat-tailed ewes with fine-fleeced rams of the Novokavkazskaya breed and Precoce (Soissonne type) in the Volgograd region (14). Due to insufficiently high indices of wool productivity, improving the shearing, wool, and meat quality and precocity was necessary. Since 1948, crossbred ewes have been crossed with rams of the Caucasian and, to a lesser extent, Grozny breeds. Today, sheep of the Volgograd breed are large, well-developed animals with pronounced meat forms and simultaneously high rates of wool and meat productivity (12). They are also quite precocious; the body weight of lambs reaches 30-35 kg before weaning, and young ewes reach up to 80% of the weight of an adult sheep by one year. Rams aged 7-9 months have a net carcass weight of 20-24 kg. The wool of Volgograd sheep is white and relatively thick. Close fleece has a staple structure and medium density. The crimp of wool fibers is quite pronounced, uniform, and somewhat stretched. Wool fineness is from 20.6 to 25.0 microns. The length of wool is 8-9 cm in ewes and 9.5-10.5 cm in rams. The pure fiber yield is 48-50%. Suint is light in color, primarily light cream. However, this breed has high-quality wool; it could be drier and more untrue. For this reason, the Volgograd sheep breed is best suited for those sheep breeders specializing in meat production. The best breeding herds have been preserved in the Volgograd region. However, the Volgograd breed can be found throughout central Russia, for example, in the Urals and the Volga region (Figure 2).

Figure 2: Volgograd sheep breed.

Kalmyk fat-tailed sheep were brought to Russia by Kalmyks in the 17th century from Mongolia and Western China during resettlement (14). As the Kalmyks moved, part of the sheep mixed with local fat-tailed sheep and formed crossbred populations, resulting in a separate Edilbay breed in Western Kazakhstan. There are some shared features of the constitution of Kalmyk and Hissar sheep. From 1930 to 1990, the accumulation of crossbreeding of fat-tailed sheep with sheep of fine-wool breeds was practiced on a large scale in Kalmykia. This led to fat-tailed Kalmyk sheep being bred only in private farmsteads and small farms. There were no pedigree farms for breeding the Kalmyk sheep breed in the republic. Since the 90s of the last centuries, an urgent need for developing early maturing mutton-fat sheep inspired an intensive revival of the Kalmyk fat-tailed sheep as a unique gene pool of aboriginal sheep breeding in the Republic of Kalmykia. The sheep have the most valuable economic and biological traits, i.e., high resistance, adaptability to a sharp temperature drop, unpretentiousness, adaptability to year-round grazing, and a combination of good quality sheepskin and meat productivity (12,17). The Kalmyk fat-tailed breed is one of the most early maturing breeds. The lamb gains body weight from 43 to 52 kg by 6-8 months of age (Figure 3).

Figure 3: Kalmyk sheep breed.

For the research, 192 individuals were selected, i.e., 42 sheep of Kalmyk fat-tailed breed (non-public joint stock company breeding plant "Kirovsky," Yashkul district, the Republic of Kalmykia), 100 sheep of the Edilbay breed (LLC "Volgograd-Edilbay," Bykovsky district, the Volgograd region) and 50 sheep of Volgograd breed (agricultural production co-operative breeding plant "Romashkovsky," Pallasovsky district, the Volgograd region). Blood sampling and animal handling were practiced according to the ethical guidelines of the Federal Research Center for Animal Husbandry named after Academy Member L.K. Ernst, the guidelines for keeping animals established in the farms, the Guide for the Care and Use of Laboratory Animals (18), and local animal care laws and policies. The current study complied with the ethical standards; the sheep owners provided written consent for using animals in this study. Blood sampling was performed by taking 5 ml of blood from the jugular vein and putting it into EDTA tubes. The tubes were stored immediately in Iceboxes, and shortly after, they were stored at −20°C until the DNA extraction step was performed.

DNA extraction

DNA was extracted from the collected blood samples using the phenol-chloroform method. DNA quality was assessed using a set of reagents: QuantiFluor dsDNA System E2670 (Promega Corp., USA) and the Quantus fluorometer (Promega Corp., USA), which measures the concentration of nucleic acids. The OD260/OD280 ratio of DNA solutions was determined by NanoDrop-2000 (Thermo Fisher Scientific, Wilmington, DE, USA).

Quality control of genotyping and preparation of genotypes

The overall genotyping quality of the obtained samples was 99.4%. The following SNP filtering parameters were selected and used in the research work; criterion of minor allele frequency (MAF). 1116 SNPs with less than 1% frequency were removed. Hardy-Weinberg equation test. 444 SNPs with a confidence threshold higher than 0.001 were removed; - test of skipping genotypes, when a single marker must be successfully genotyped in at least 90% of the animals in the sample. 824 SNPs genotyped in less than 10% of the population were excluded. Processing raw information resulted in 47882 markers on 26 chromosomes left for further work (Figure 4).

Figure 4: Distribution of SNP markers on sheep chromosomes.

Genotyping and data analysis

Genotyping was performed using the OvineSNP50 DNA Analysis Kit (Illumina Inc., USA). Each genotype contained 50266 alleles of single nucleotide polymorphisms (from now on referred to as SNP). Genotyping quality control, inbreeding analysis, and IBS kinship matrix construction were conducted using Plink 1.9 software (19). Model clustering was carried out using the Admixture program (20). The association analysis applied to the EMMAX program. The analysis results were graphically visualized using the R software environment (21).

Identification of SNPs associated with body weight

Based on the body weight data, an association analysis was performed to search for genes and genome regions associated with the body weight of sheep. The analysis was conducted in the EMMAX program using a mixed model equation that has the form of y = µ + Br + Xb + e. Where y is the body weight indicator, kg; Br is the breed of sheep; Xb is the SNP marker effect; and e is the residual unknown. The p-value was calculated for each marker.

Results

Assessment of molecular inbreeding

Regarding individual indicators of molecular inbreeding, the average level of inbreeding was 2.72% for the Volgograd breed, 1.96% for the Kalmyk breed, and -0.0022 for the Edilbay breed. A negative indicator of genomic inbreeding indicated a high molecular heterogeneity of the herd under study (Figure 5).

Figure 5: Genomic inbreeding coefficient of the studied breeds. (A) Volgograd; (B) Kalmyk; and (C) Edilbay breeds.

Assessment of genetic diversity in genotyped animals

Genetic diversity was assessed to detect genetic differences between breeds and intrabreed groups. According to the data obtained, two a priori known breeds-Volgograd (red dots) and Kalmyk (green and blue dots)-were distinguished. The location of the breeds on the graph indicated the absence of a genetic relationship between the two breeds represented (Figure 6).

The Volgograd breed of sheep was divided into three subbreed groups (clusters). Central cluster (A) was the most numerous modern sheep-bred group se. Clusters B and C were separate along the y-axis, indicating a distinctive genetic structure compared to the central cluster. This division may be due to the historical admixture of blood of Grozny and French Soissons breeds. A detailed study of the sheep pedigree in the clusters allowed for the correct planning of breeding work with these animals. The Kalmyk sheep had a single structure and were not distinguished into additional subclusters. Two individuals marked as cluster D on the graph had increased genomic inbreeding coefficients. For this reason, they were located remotely from the main array. It was advisable to exclude (cull) these sheep from the breeding work in the herd.

The Edilbay sheep were divided into three clusters (Figure 7): (A) and (C) were separately identified according to the distribution of sheep along the abscissa; at the same time, in cluster (C), two separate subclusters (C1) and (C2) were additionally distinguished; along the ordinate axis, a group of animals named as cluster (B) is distinguished from cluster (A) was identified. A detailed study of the pedigree of sheep from different groups made it possible to understand the reason for the genetic differences of individuals and correctly plan the breeding work. The result of multidimensional scaling of all genotyped individuals is shown in figure 8.

According to the abscissa axis, the genetic distance between the Kalmyk and Edilbay sheep breeds was significantly lower than between the Volgograd and Kalmyk or Edilbay breeds. The graphical analysis indicated a high relationship between the Kalmyk and Edilbay sheep breeds (Figure 8). Based on the results shown earlier (Figure 6), the Volgograd breed was divided into three breed groups along the y-axis. The information obtained was of practical importance for selecting and breeding work on these farm animals.

Figure 6: Graphical display of multidimensional scaling based on genome-wide data of the studied breeds: red dots are the Volgograd breed; green and blue dots are the Kalmyk breed. (A), (B) and (C) are the subclusters of the Volgograd breed. (D) are two individuals with an increased inbreeding coefficient (Kalmyk breed cluster).

Figure 7: Graphical display of multivariate scaling based on genome-wide data of the Edilbay breed.

Figure 8: Graphical display of multivariate scaling based on genome-wide data of three breeds: red is the Volgograd breed; green is the Kalmyk breed; and blue is the Edilbay breed.

Genetic features of Volgograd and Kalmyk breeds analyzed by the method of model clustering

A cross-validation (CV) procedure was performed at the first stage to determine the a priori population number (K) in the analyzed sample. According to the procedure's results (Figure 9), the minimum CV value was observed at K=2. Thus, the algorithm assumed the presence of two contrasting populations (breeds) in the presented sample. Figure 10 shows a graphical division of sheep into two populations according to the procedure of the model clustering method. At K=2, the number of populations was equivalent to the number of analyzed breeds. According to the multidimensional scaling results, the samples were divided into four subpopulations, i.e., one cluster of the Kalmyk breed and 3 clusters of the Volgograd breed. For this reason, the current sample was analyzed with the criterion K=4 (Figure 11). The obtained MDS graph clearly showed that separated clusters B and C were monochrome colored in green and blue, respectively. Some Volgograd sheep (cluster A) had spots of these colors in their structure. This made it possible to judge the potential belonging of the samples in clusters B and C to the original breeds or types. Cluster A samples were a mixture of a third breed (red) with sheep in clusters B and C.

Figure 9: Results of the cross-validation (CV) procedure to determine the a priori population number (K) in the analyzed sample.

Figure 10: The distribution graph between the two races.

Figure 11: The distribution graph between the two races considering three clusters of Volgograd sheep breed.

The genetic features of the Edilbay breed were analyzed using model clustering

In the first stage, a CV procedure was performed to determine the population number (K) in the analyzed sample. According to the procedure's results (Figure 12), the minimum CV value was found at K=1. In the presented sample, the algorithm did not allow for the identification of differences between the clusters or breed types compared with interbreed analysis.

Figure 13 graphically shows the studied sheep’s genetic structure identified by the model clustering method. According to the P.C. analysis results, K=5 corresponded to the number of clusters and was chosen as the a priori value. Most of the samples had a mixed genetic architecture in different proportions. Separate 18 individuals were characterized by a unique, stable monomorphic genetic structure (painted in only one color). A similar staining color of the monomorphic individuals indicated an identical genetic structure of these samples. According to the calculation results, 18 sheep formed the following monomorphic groups: red, blue, green, light blue, and yellow. The unique structure of these sheep should be considered in the breeding work. The fact that the rest of the sheep were stained with the same colors but in different proportions enabled judging the 18 monomorphic samples belonging to the original groups.

Figure 12: Results of the CV procedure to determine the a priori population number (K) in the analyzed sample (К).

Figure 13: The distribution graph of Edilbay breed individuals.

Complex analysis of genetic features of three breeds by the model clustering method

The lowest CV coefficient value was obtained at K=3 (Figure 14), and an addition index decrease was observed at K =5. The presence of three subgroups in the Volgograd sheep breed substantiated this value. Figures 15 and 16 show a graph that allowed the races to be divided into 3 and 5 groups, respectively. In both cases, large sections of the genome of the Kalmyk and Edilbay races were stained in similar colors. Moreover, the Volgograd breed was genetically isolated.

Figure 14. Results of the CV procedure for determining the a priori number of populations in a sample consisting of 192 head of the Volgograd, Kalmyk, and Edilbay breeds of sheep.

Figure 15: The distribution graph between the three races.

Figure 16: The distribution graph between three races considering three clusters of Volgograd sheep breed.

Identification of SNPs associated with body weight.

Figure 17 shows a distribution graph of confidence values (p-value) for markers distributed on chromosomes. The highest relationship was found for the marker OAR10_26672645.1 (p=4*10^-8) located on chromosome 10. In the chromosome region 50.26-51.51Mb (Figure 18), there are no genes that have been described had a potential relationship with this trait. Another marker that passed a significant confidence threshold was OAR1_18293636.1 with p=7.3*10^-7. Two previously described genes, FARP2 and HDLBP, were located on chromosome 1 in the region of 32.69 - 54.87 Mb (Figure 19). Of most significant interest was the HDLBP (High Density Lipoprotein Binding Protein) gene, which is responsible for synthesizing protein that binds high-density lipoprotein and the level of cholesterol in the cell.

Figure 17: The result of a genome-wide association analysis with the trait “body weight” in the sheep population.

Figure 18: The chromosome 10 region fragment that shows the highest significance value in the association analysis.

Figure 19: The chromosome 1 region fragment that shows the highest significance value in the association analysis.

Discussion

Evolution studies to understand the genetic potential of animals and breeds within animal production are of great value, as they generate good results that can be applied not only to improving the species but also to the economic improvement of the activity performed. As known, a relatively close relationship between the Edilbay sheep breed and the Kalmyk fat-tailed sheep was confirmed by only the microsatellite analysis of sheep DNA. Minimum genetic distances were also established pairwise in Grozny-Edilbay, Grozny-Karakul, and Edilbay-Karakul breeds, i.e. 0.1364 and 0.0851; 0.1620 and 0.1208; 0.1875 and 0.1192, respectively (22). Besides, in other regions of Russia, 25 local populations of sheep breeds, including Volgograd (fine-fleeced), Edilbay, and Kalmyk (coarse-wool), were studied. The research was conducted using 11 microsatellite loci (OarCP49, INRA063, HSC, OarAE129, MAF214, OarFCB11, INRA005, SPS113, INRA23, MAF65 and McM527) (23). Our study-based genome-wide analysis also noted a high relationship between the Kalmyk and Edilbay sheep breeds in the studying regions. At the same time, the Volgograd breed was genetically isolated. Additionally, for the first time, we found that the Volgograd breed of sheep was divided into three separate subbreed clusters because of the historically targeted admixture of blood of other, more productive sheep breeds. Pooling obtained genetic characteristics of the major Russian local sheep breeds that are moderately diverse and have a strong population structure with a worldwide genotyping set, which gave deeper insight into the history and origin of the Russian sheep populations. For example, the Russian fat-tailed breeds shared co-ancestry with sheep from China and Southwestern Asia (Iran) (14). Also, the genetic diversity, structure, and shared ancestry between sheep breeds found in the Kingdom of Saudi Arabia (KSA) were investigated by analyzing genetic variation in microsatellite markers (9). The sheep breeds of the KSA revealed high genetic diversity even though they are reared in different geographic regions that are far apart and have different features. In Kazakhstan, 75 individuals from five sheep populations were investigated based on 12 STR (short tandem repeat, also known as microsatellite) markers to study their genetic structure and phylogenetic relationship based on genetic distances (24). As a result, the highest genetic distance was observed between Kazakh Arkhar-Merino and Edilbay-1, whereas the genetic distance between Edilbay-1 and Edilbay-2 is the smallest using Nei's standard genetic distance. The Edilbay-1 sheep breed possesses the largest genetic diversity among the five studied populations. According to the hypothesis, the Edilbay breed animals could also be presented by the original and new breed types, marked on the graph in green and red. However, according to the results obtained, dividing samples into two clusters according to the breed type was impossible because all values were in admixture throughout the graphic field. For the first time, a medium-scale SNP genotyping of five local sheep breeds in Kazakhstan was applied to investigate their population structure and relation to global sheep diversity (25). Principal component analysis and model-based structure analysis of general population markers revealed two breed groups. The first group included Akzhayik and Kazakh fine-wool sheep; the second group included Edilbay, Saryarka, and Kazakh semi-coarse wool sheep. The heterogeneity of different Akzhayik and Kazakh Semi-coarse wool sheep populations was observed. A neighbor-joining tree comparing Kazakh sheep data with the dataset generated by the Sheep HapMap project supported a close relationship between Kazakh sheep varieties and ancient domestic sheep ancestors. A comprehensive study of genetic characteristics of sheep breeds in Great Britain and their importance for the global gene pool due to tools and applications for assessing, breeding, and producing livestock, including, above all, genetic profiling of various breeds and searching for quantitative trait loci (QTLs) and candidate genes in farm animals was conducted (3). British sheep breeds are common in other parts of the world and have been exported from the U.K. to other countries to create new breeds and improve existing ones, with Australian, American, and European breeds being prime examples. One can mention imports to the former USSR, when British breeds served as the basis for new breeds, such as the Gorky and Russian long-wool breeds. In particular, British sheep meat breeds noticeably affected the development of local breeds in Russia and the countries of the former USSR in the last century. The article discusses the uniqueness of British breeds in terms of phenotypic traits and adaptation to local conditions. Also, a study to identify the genetic characteristics of three sheep breeds from Azerbaijan, namely, Bozakh (n = 19), Karabakh (n = 16), and Mazekh (n = 19), was conducted (10). The studied Azerbaijani sheep breeds were established to be genetically close but carry specific genomic components, making it possible to consider them valuable genetic resources. So, in the evolution, all breeds formed their clusters, while Bozakh was obtained from the Karabakh and Mazekh breeds. Pairwise, FST distances were quite low, with the highest value of 0.017 between the Karabakh and Mazekh breeds and the lowest value of 0.006 between the Mazekh and Bozakh breeds. The global fixation index (FST) was 0.01, indicating that 99% of the genetic variation was due to differences within breeds. The expected heterozygosity was the same in Bozakh and Mazekh (0.363±0.000 and 0.360±0.000, respectively) and somewhat lower in the Karabakh breed (0.353±0.000). Neither excess nor deficiency of heterozygotes was observed in all the breeds studied. The inbreeding coefficient (Fis) ranged from 0.001 in Karabakh to 0.011 in Mazekh sheep. It was found that the Gotland sheep breed (Sweden) is a particular breed but also has shared ancestral genomic components with Gute (~50%), Karakul (~30%), Romanov (~20%), and Fjällnäs (~10%) breeds (11). Thus, it is essential to investigate the genetic features of native and foreign breeds of sheep by different modern methods.

In the previous studies based on SNPs in significant association with three milk production traits (milk yield, fat yield, and protein yield) in a crossbred dairy sheep population by three statistical techniques, including GWAS, were identified (6). Later, the association between the estimated breeding values of the milk yield with 50 K SNP chip results in Awassi sheep was investigated (7). The results of a genome-wide association analysis in a resource sheep population (Ovis aries) originated from backcrosses of (Romanov ½ Katahdin) ½ Romanov were also presented (13). Their body weights were recorded in age dynamics, and SNP profiles were obtained using a high-density DNA chip. The GWAS-analysis of the body weight of backcrosses (Romanov ½ Katahdin) ½ Romanov revealed 38 SNPs significantly associated with the body weight (p<0.00001) on OAR1, OAR2, OAR3, OAR4, OAR6, OAR9, OAR10, OAR11, OAR13, OAR15 and OAR19, as well as functional candidate genes affecting skeletal muscle growth, bone skeleton formation and lipid and carbohydrate metabolism. In addition, the authors showed age-related changes in the composition of significant SNPs. Our current studies have also shown the highest relationship between the body weight with the markers OAR1_18293636.1 (p=7.3*10-7) and OAR10_26672645.1 (p=4*10-8).

Conclusion

The performed study contributes to a deeper understanding of the mechanisms and patterns of the genetic structure of sheep populations developed in specific agroecological conditions in the south of Russia in the evolution, consolidation of economically important phenotypic traits in the ontogenesis and an increase in local breed's potential fulfilled, which is especially important when solving the problem of intensifying the sheep breeding industry. A detailed study of the sheep pedigree in the selected subclusters will make it possible to understand the reason for the genetic diversity of individuals and correctly plan the breeding work with these animals.

Acknowledgment

The authors are grateful to the Russian Science Foundation for financial support of this study (Project 22-16-00041 NIIMMP), and owners and employees of farms where the investigation was conducted are also to.

Conflict of interest

There is no conflict of interest.