Genome-Wide Association Study (GWAS) For Left Displaced Abomasum in Highly Productive Russian Holstein Cattle

1. Introduction

Left displaced abomasum (LDA) is a multifactorial disease of cattle that occurs mainly during the transition postpartum period [1,2]. In the left form of displacement, the abomasum occupies a caudodorsal position, located between the rumen and the left abdominal wall. The pathology is characterized by a decrease in milk production and an increased risk of culling. The most effective way to treat this pathology is surgery, which significantly increases the economic costs of milk production. The median incidence of left displaced abomasum among Holstein cattle based on 12 studies is 2.71% [3], which makes this pathology one of the important threats to the economic stability of dairy livestock farming.

Numerous studies have established that left displacement of the abomasum is recorded often (75% of cases) than right displaced abomasum (RDA) [4,5]. Right displacement is more likely to be persistent, in contrast to left-sided form, which significantly increases the risk of right pathology. It has also been established that right form is more acute than left one and requires prompt surgical treatment. However, due to its frequent occurrence, the economic damage from LDA significantly exceeds the losses from RDA.

Despite the widespread occurrence of LDA, its etiology remains not fully understood. There are several hypotheses explaining LDA etiology. It is believed that one of the most important prerequisites for displacement of the abomasum is atony arising as a result of inadequate feeding that does not correspond to the physiological status and leading to the accumulation of gas and liquid in the organ cavity, since this symptom occurs in more than 80% of all cases [4,6,7]. According to the hypothesis, gas appears in the abomasum body and, with an increase in its volume, shifts abomasum cranially. The flow of gas toward the pyloric region causes displacement of the abomasum to the cranial dorsal side and then its movement to the left side of the abdominal wall [8]. Occurrence of left displaced abomasum was also associated with the elliptical shape of the abdominal cavity section [9] and negative energy balance [10]. It has also been demonstrated that an increase in BHB and AST levels in the first two weeks after birth is a predictor for the occurrence of the described pathology [11,12,13], as well as low concentrations of liver-secreted IGF1 together with increased non-esterified fatty acids (NEFA) and beta-hydroxybutyrate (BHB) [13]. Cows with left displaced abomasum often have comorbidities such as ketosis, mastitis, metritis, endometritis and fatty liver [4,14].

Heritability for LDA has been estimated in Holstein cattle to be between 0.1 and 0.31 [2,15,16]. Early attempts have already been made to determine the genetic determinants for the disease. Recent genome-wide association studies (GWAS) have identified several genes and single nucleotide polymorphisms (SNP), associated with left displaced abomasum [2,17,18]. However, clarifying the position of quantitative trait loci (QTL) in individual populations and identifying the genes underlying this complex disease can contribute to both a deeper understanding of the pathogenesis of LDA and the development of effective approaches to the treatment and prevention of this pathology.

Another effective tool for genomic analysis is the assessment of individual autozygosity of animals, or extended homozygous nucleotide fragments - runs of homozygosity (ROH), indicating the inbreeding. The consolidation of ROH loci with a low frequency of recombination in generations of animals may indicate the formation of gene clusters under selection pressure, mostly due to positive selection for productivity traits, as well as adaptation to changing environmental conditions. ROH allow us to more effectively assess the degree of genomic inbreeding both in a population and individually and also to evaluate the conservation of certain loci. A high degree of inbreeding in the studied animals can influence the results of genome-wide association study, increasing the likelihood of detecting false positive associations. Moreover, the location of candidate genes in ROH may indicate their significant influence on the breeding quality of animals. Modern industrial dairy farming technologies are based on the exclusive use of assisted reproductive technologies, such as artificial insemination of cows. For cryopreservation of genetic material (sperm), sires with the highest breeding value and high rates of heritability of productive qualities are selected. This approach leads to an increase in the inbreeding degree in herds and an increase in the population genetic homogeneity. As a result of the use of such a strategy for improving the genetic potential of dairy cattle, the frequency of predisposition hereditary factors for multifactorial diseases increases. Loci containing disease susceptibility factors often become homozygous, which increases the number of disease cases in herds.

The aim of the study is to find genetic associations characterizing the genomic variability of the predisposition to LDA in Holstein cattle of the Leningrad region of the Russian Federation.

2. Materials and Methods

The collection of biological material for the study was carried out in three farms in the Leningrad region. The study analyzed two groups of animals, divided according to the principle of “case” and “control” (healthy animals). Peripheral blood samples were taken from a mixed-age herd of 360 highly productive Holstein cows, from which 175 animals belonged to the control group and 185 cows belonged to the case group that consisted of animals with left displaced abomasum diagnosed by clinical examination. Sampling was carried out in specialized vacutainer blood collection tubes with EDTA as an anticoagulant then stored in frozen condition at -200C until the DNA extraction procedure. DNA was isolated from samples using the phenol-chloroform method according to a standard proteolytic treatment protocol using proteinase K. Quality control of isolated DNA was carried out using a NanoDrop 2000 spectrophotometer (Thermo Fisher, USA).

For the genotyping procedure the GeneSeek® Genomic Profiler™ Bovine 100K DNA chip (Neogen, USA) with a coverage density of 95256 SNPs was used. Principal component analysis (PCA) was performed to assess population stratification and determine the number of genetic covariates. Quality control and filtering of genotyping data were carried out using the PLINK 1.90 software package [19]. When filtering the bioarray data, the following value thresholds were used: the minimum allele frequency was set at a threshold of <0.5, the p-value threshold for the Hardy-Weinberg distribution was <1e-6, the number of missing variants per sample was <0.05, and the proportion of missing variants per marker was also <0.05.

ROH analysis was implanted in detectRUNS package in R. Parameters for ROH search were determined according to the following criteria [20]: a sliding window size was 25 SNP, the proportion of homozygous overlapping windows was 0.05, a minimum number of 35 homozygous SNPs included in a window, a maximum gap between consecutive SNPs of 0.1 Mb, a maximum number of 1 heterozygous SNP per window. Found ROH islands then were divided into three classes based on its size: short (0.5-1Mb), medium (1-4Mb) and large (>4Mb) [20]. To evaluate the inbreeding coefficient, we count proportion of the genome covered by runs of homozygosity by dividing total length of ROH by the length of genome (FROH) [21] according to the following formula:

$${{\mathrm{F}}_{\mathrm{R}\mathrm{O}\mathrm{H}}}_{\mathrm{i}}={\displaystyle \frac{{\sum }_{k=1}^n{{\mathrm{L}}_{\mathrm{R}\mathrm{O}\mathrm{H}}}_i}{{\mathrm{L}}_{\mathrm{g}}}}$$

(1)

where:

• i is the number of the animal_i,
• n stands for total number of ROH found in animal_i,
• L_ROHi represents total ROH length for animal_i,
• L_g identifies genome length of the individual.

We defined a common ROH island based on overlapping homozygous regions with a frequency greater than 0.50 for all animals and also considered the suggestive frequency threshold greater than 0.40 to include more candidate genes and regions. Сandidate genes located within ROH islands were identified using cattle genome assembly ARS-UCD1.3 (https://www.ncbi.nlm.nih.gov/datasets/genome/GCA_002263795.3/). To find genomic loci that may contain genes whose polymorphisms are strongly associated with the development of LDA, we also looked at the distribution of ROH islands in the case and control groups separately. The statistical significance of the difference in the incidence of ROH islands in healthy animals and cows with displacement was assessed using the Pearson chi-square test.

To ensure accuracy of genome-wide association study we used three statistical models independently: General Linear Model (GLM), Mixed Linear Model (MLM), and Fixed and Random Model Circulating Probability Unification (FarmCPU), implemented in the GAPIT Version 3 [22] and rMVP [23] packages in R.

False discovery rates (FDR) were calculated as described by Benjamini and Hochberg [24]:

$$\mathrm{F}\mathrm{D}\mathrm{R}={\displaystyle \frac{\mathrm{m}\times \mathrm{P}\mathrm{r}\left(>0\right)}{\mathrm{R}}}$$

(2)

where:

• m is the number of tests,
• R is the number of significant tests,
• Pr(> 0) is equal to the significance threshold.

The results were visualized using the qqman (v.0.1.9) [25] and qqplot2 (v.3.4.4.) packages in the RStudio programming environment (v. 2023.9.1.494) using the Manhattan plot.

Annotation of the identified associations was carried out using cattle genome assembly ARS-UCD1.3, visualized in the genome browser GDV (Genome Data Viewer). To search for QTL associated with Economically relevant traits the CattleQTLdb database was used. For functional annotation of genes located within ±0.2 Mb of suggestive SNPs the National Center for Biotechnology Information (NCBI) Gene database was used, as well as the open UniProt protein sequence database and scientific literature. The impact of a substitution on a protein sequence was determined using the VEP (Variant Effect Predictor) software [26].

3. Results

In our study 85,605 polymorphic genetic variants obtained from genotyping by GeneSeek® Genomic Profiler™ Bovine 100K DNA array and 325 samples passed quality control. Of the filtered samples 172 were taken from animals of experimental group and 153 of the control group.

We evaluated inbreeding coefficient using F_ROH method. The F_ROH ranged from 0.01 to 0.19 for individual genotypes and from 0.09 to 0.17 for mean F_ROH value for each bovine chromosome. We found that the distribution of homozygosity has similar value across autosomes. However, BTA 9 and BTA 15 are distinguished by a noticeably lower content of runs of homozygosity, compared to other chromosomes. On the contrary, BTA 20 in animals from our sample on average contains more homozygous regions, which requires further consideration. The distribution patterns of genomic inbreeding coefficients among the studied animals and bovine chromosomes are shown on Figure 1.

Analysis of ROH islands occurrence frequency in the studied population of dairy cattle allowed us to identify a number of conserved regions in the genome, which, as we assume, are under selection pressure. As a next step, we calculated the frequency of ROH and identified ROH islands for each chromosome in our population. First, we identified regions with an occurrence frequency greater than 0.5. Two such ROH island was discovered on BTA14 and two were located on BTX (Table 1).

Also, 20 homozygous regions were found, including at least 5 SNPs, which were found in more than 40% and less than 50% of the all studied population. Suggestive ROH were mapped to BTA5, BTA7, BTA13, BTA14, as well as BTA20, BTA21, BTA22, BTA26 and BTX. Suggestive homozygous regions on BTX had the longest extend (6.4 Mb, 5.5 Mb and 3.9 Mb) and were present in the greatest number (n = 3) which can be explained by the characteristics of sex chromosomes recombination.

The distribution of ROH islands in the control and case groups of animals was slightly different. In particular, four ROH, two of which were located on BTA13, one on BTA21 and one on BTX were more frequently represented in the genomes of animals in the control group. It is worth noting that no selection signatures were found that were significantly more common in the control group than in the experimental group. Information about extent of the loci and their frequency in the both groups is shown in Table 2.

According to the CattleQTLdb first loci on BTA13 contains QTL associated with body conformation traits such as stature, dairy form, rear leg placement, body depth, rump width, foot angle and few traits related to productivity and susceptibility, namely milk yield and milk fat yield and M. paratuberculosis susceptibility. The second ROH on the given chromosome, which has the largest frequency range among the two groups, is also contains QTL for a similar set of traits, specifically body depth, rump width, stature, dairy form, udder cleft and several others. Moreover, the occurrence of this region in the genomes of animals of the two groups was significantly different, which indicates the role of body shape in the development of the disease.

The region of homozygosity on BTA21 does not contain the described QTL but in humans, loss of homologous genomic region is connected with the Prader-Willi syndrome. ROH island on BTX is associated with semen quality more specifically percentage normal sperm and general reproductive traits such as age at puberty and scrotal circumference.

Regarding the GWAS analysis, the choice of three models was based on their unique advantages in addressing various complexities of genetic analysis. GLM is an effective tool for analyzing simple population structures, whereas MLM provides a more robust framework for dealing with complex population structures, potentially reducing false positive associations [27]. FarmCPU takes a novel approach implementing both fixed and random effects, allowing to simultaneously handle effects of cryptic relatedness and population stratification thus make it possible to more accurate map quantitative trait loci [28]. Our comparative analysis revealed subtle differences in the results of all three models. Thus, we focused on the MLM and FarmCPU models as this could minimize the chance of finding false associations while increasing the likelihood of detecting true associations, resulting in a more robust and complete set of marker trait associations.

Figure 2 shows the results of a genome-wide association study with left displaced abomasum in the studied group of Holstein cows. Fourteen suggestive SNPs (p < 1 × 10^-4) were detected. The found polymorphisms were identified on BTA2, BTA3, BTA6, BTA9, as well as BTA16, BTA23, BTA24, BTA25, BTA27 and BTX.

Figure 1. Manhattan plot of significant SNPs for left displaced abomasum in Russian Holstein cattle identified in FarmCPU, GLM, and MLM of GWAS analysis. The lower horizontal line represents the suggestive significance threshold (p < 1 × 10⁻⁴). The upper line represents genome-wide significance threshold (p < 1 × 10⁻⁴). Colored dots indicate an associated SNP.

Figure 1. Manhattan plot of significant SNPs for left displaced abomasum in Russian Holstein cattle identified in FarmCPU, GLM, and MLM of GWAS analysis. The lower horizontal line represents the suggestive significance threshold (p < 1 × 10⁻⁴). The upper line represents genome-wide significance threshold (p < 1 × 10⁻⁴). Colored dots indicate an associated SNP.

The largest number of suggestive SNPs was located on BTA25 (n = 4); all polymorphisms were located within locus positioned in range from 17.7 to 27.1 Mb. Two nucleotide polymorphisms were located on BTA16, the remaining chromosomes each contained one nucleotide substitution associated with the development of pathology.

The structural annotation revealed 14 genes that contained detected SNPs or were close to them. Functional annotation of polymorphisms showed the presence of substitutions mainly in the non-coding parts of the genome (such as intergenic regions, introns). Of the 14 polymorphisms detected, of which 6 were in introns of following genes: ABCB11, SRP72, RGS18, GSG1L, FBXL19 and PRKCB. Remaining polymorphisms were found in intergenic regions. Substitutions located in exons have not passed the threshold of statistical significance, were synonymous and did not affect the functions of proteins encoded by the genes in which they were located. More detailed information about the associated polymorphisms is in Table 3.

Two SNPs were located in QTLs previously described for the trait. Thus, the SNP located on BTA23 was within the boundaries of the locus with coordinates 29.2-47.0 Mb [29]. Also, in previous studies on the Chinese population of Holstein cattle, one of the regions on chromosome 25 (24.8–26.8 Mb) was already noted [Huang H].

GO enrichment analysis showed that candidate genes were involved in a wide range of biological processes, including protein ubiquitination (ABCB11, FBXL19), fatty acid metabolism (ABCB11, PNPLA4), the formation of organs and tissues as homeotic genes (SOX4), transcription regulation, regulation of signal transduction and others. A list of GO (Gene Ontology) terms characterizing biological processes is displayed in the table below.

However, the biological role of certain genes has not been established due to their low representation in the scientific articles.

4. Discussion

Genome-wide association study for left displaced abomasum made it possible to identify associated polymorphisms localized close to or within genes, of which such genes asABCB11, SRP72, RGS18, SOX4, GSG1L, FBXL19, PNPLA4 are most promising. Our results suggest that susceptibility to LDA may be formed by the contribution of multiple genes with small effect sizes, consistent with moderate estimates of heritability for the trait.

Past study of Ricken et al. have found genetic correlations between LDA and milk fat, milk protein and milk yield [30]. Our results partly confirm previous findings, as some of the genes described in our study located in QTL for traits such as milk fat content (ABCB11) [31] and glycerophosphocholine content (GSGL1) [32].

Runs of homozygosity analysis shown that the studied population is characterized by a moderate degree of genomic inbreeding. For most animals this value ranged from 0.1 to 0.15. Four ROH islands were discovered that were present in the genomes of most animals of both groups, all of them located on two chromosomes - BTA14 and BTX. Most likely, the reason for the appearance of such regions in the genome is selection pressure aimed at improving the quality of livestock productivity. It should be noted that some of the discovered regions, as well as the candidate genes located in them, have already been described for other herds in the scientific articles. As for one of homozygous region located on BTA14, genes located in the locus have already been associated with such phenotypic traits as serum prolactin level [33], growth and feed efficiency [34], carcass weight (CWT) and eye muscle area (EMA) in Korean Hanwoo cattle [35], carcass trait [36], fat deposition in Brazilian Nellore cattle but productive traits, especially important for dairy cattle, have not yet been associated with this locus.

Some loci were more widely represented in the group of animals with LDA than in healthy individuals. Thus, two loci on BTA13, both containing QTL for body conformation traits had the greatest difference in prevalence (10.11% and 15.28%, respectively) between two groups. Moreover, the second ROH was statistically significantly different in its occurrence in the case and control groups and was more widely represented in animals with LDА. That may be due to them containing selection signatures referring to the body conformation of Holstein cows. As previously mentioned, the occurrence of LDA was associated with elliptical shape of the abdominal cavity section [9] which is typical for this breed of cattle. The shape of the abdominal cavity may affect the location of the internal organs, causing the abomasum to be pushed out of the right side of the body when it is inflated. Although body measurements were not taken in this study, we can hypothesize that the shape and volume of the abdominal cavity, according to genomic data, may be a predisposing factor for left displaced abomasum in Holstein cattle.

Another ROH more common in animals with LDA was located on BTA21. In our sample, healthy cows were 7.59% less likely to become its carriers. The locus contains five copies of genes coding small nucleolar RNASNORD116 and one copy of SNORD109A. In humans, this genes strongly associated with Prader-Willi Syndrome and Angelman Syndrome, complex multisystem genetic disorders with developmental delay, behavioral uniqueness and speech impairment [37,38]. Prader-Willi Syndrome is also characterized by life-threatening obesity due to uncontrolled appetite and muscle weakness [38]. Despite the fact that QTL has not yet been identified for productive traits in this genetic region, it may play a significant role in the formation of the nervous system and appetite control both in humans and cattle according to the principle of orthology which may be connected with occurrence of the studied disease.

ABCB11 encodes the bile acid export pump (BSEP), the main transporter of bile acids from hepatocytes to the biliary system [39]. The connection between liver diseases and impaired motility of the digestive system in cattle has been demonstrated in a number of publications [2,4,40,41]. Thus, in cows with pathologies of this organ, a decrease in the strength and duration of rumen contraction was noted [41]. Studies have reported that cases of LDA in animals with liver lesions are more common than in the general population [2]. When analyzing the biochemical profile of animals, it was found that in cases of LDA, markers of hepatocyte death significantly increase [4]. And, vice versa, an increase in some of hepatocyte death markers predicts the onset of displacement [11,12], however, to explain this connection requires further research. A previous GWAS study has already suggested a link between bile acid metabolism and LDA [2].

The SRP72 gene encodes a polypeptide that is part of a protein complex that promotes post-translational modification of proteins in the endoplasmic reticulum [43]. In the signal recognition particle (SRP) SRP68 and SRP72 proteins form a heterodimer that has been difficult to investigate and thus many processes in which the previously mentioned protein is involved still remain not fully understood. It is known that polymorphisms in SRP72 are associated with one of the hereditary forms of aplastic anemia [44]. Transcriptomic studies have demonstrated that the expression of SRP72 gene may be a marker of hepatotoxicity [45]. The mechanism of involvement of this protein in the development of left displacement of the abomasum is also not completely clear. SRP72 may carry out post-translational modifications of proteins, whose direct or indirect interaction with certain substrates will lead to a decrease in the force of smooth muscles contraction which may affect on peristatic of the gastrointestinal tract and cause the occurrence of LDA. Or the effect may be more indirect, for example, it can be caused by changes in the passage of certain biochemical pathways.

SOX4 encodes a transcription factor involved in the formation of several organ systems during embryogenesis [46]. The possible role of this gene in a hereditary predisposition to abomasum displacement is ambiguous since SOX4 is involved both in the formation of the nervous system and has a pleiotropic effect on the formation of the bile ducts of the liver [47].

The GSG1L protein is a subunit of the AMPA receptor that regulates the transmission of nerve impulses at synapses [48]. FBXL19 is involved in the ubiquitinylation of proteins and also regulates the activity of a number of homeotic genes by attaching CDK8 to their promoters [49], thereby participating in the formation of the musculoskeletal system and internal organs disturbances in this process may be a predisposition factor for LDA. The protein is also involved in IL-33 signaling pathway [50]. The role of these proteins in the development of left displaced abomasum requires further study.

Thus, we mapped loci that may be associated with LDA in Holstein cattle of the Leningrad region. Found loci were identified on BTA2, BTA3, BTA6, BTA9, as well as BTA16, BTA23, BTA24, BTA25, BTA27 and BTX. Most of the polymorphisms were not detected in previous studies, which may indicate the diversity of the genetic architecture in populations distributed in different geographical areas.