Transcription Factor Binding Site Polymorphism in the Motilin Gene Associated with Left-Sided Displacement of the Abomasum in German Holstein Cattle

Left-sided displacement of the abomasum (LDA) is a common disease in many dairy cattle breeds. A genome-wide screen for QTL for LDA in German Holstein (GH) cows indicated motilin (MLN) as a candidate gene on bovine chromosome 23. Genomic DNA sequence analysis of MLN revealed a total of 32 polymorphisms. All informative polymorphisms used for association analyses in a random sample of 1,136 GH cows confirmed MLN as a candidate for LDA. A single nucleotide polymorphism (FN298674:g.90T>C) located within the first non-coding exon of bovine MLN affects a NKX2-5 transcription factor binding site and showed significant associations (ORallele = 0.64; −log10Pallele = 6.8, −log10Pgenotype = 7.0) with LDA. An expression study gave evidence of a significantly decreased MLN expression in cows carrying the mutant allele (C). In individuals heterozygous or homozygous for the mutation, MLN expression was decreased by 89% relative to the wildtype. FN298674:g.90T>C may therefore play a role in bovine LDA via the motility of the abomasum. This MLN SNP appears useful to reduce the incidence of LDA in German Holstein cattle and provides a first step towards a deeper understanding of the genetics of LDA.


Introduction
Left-sided displacement of the abomasum (LDA) is a common dairy cattle disease with prevalences of about 2% in the German Holstein (GH) population [1,2]. In the course of LDA the abomasum starts bloating and displaces from the dexter abdominal floor to the sinistral abdominal wall. LDA has to be treated surgically by fixing the abomasum into its correct position. However, even if the treatment succeeded without any incidents, a significantly reduced milk performance as well as conception problems increase the culling risk [1][2][3]. On the average, about one half of all cows affected by LDA were culled within the first year after surgery [1]. Thus, LDA increases veterinary costs for dairy farms and decreases life expectancy for dairy cows. Several QTL for LDA have already been identified in GH cows [4], among them a QTL proximal on bovine chromosome (BTA) 23. Since LDA is usually preceded by a decreased motility of the abomasum, impaired abomasal emptying, and malfunctions at the level of the intrinsic nervous system combined with impaired cholinergic muscle responses [3], the motilin (MLN) gene located proximally on BTA23 was an appropriate candidate gene. MLN encodes a small peptide hormone regulating interdigestive gastrointestinal contraction [5]. In human, MLN is mainly expressed in the stomach (UniGene, EST profile, HS.2813), which is the analog to the bovine abomasum. A reduced expression of MLN might lead to an insufficient gastrointestinal motility. This is also supported by results in human pharmacology, where medications based on MLN-agonists can remedy conditions like the diabetic gastroparesis [6,7]. Also in cows undergoing surgical correction of LDA, the motilin agonist erythromycin increases the abomasal emptying rate by binding to motilin receptors [8].
LDA is a multifactorial disease, in which environmental effects play a role such as twin births [1,9], housing system [1], and feeding factors [10]. However, the heritability of LDA reported in various studies ranged between 20% and 50% [1,2] and thus is much higher than in any other common cattle disease.
The objective of this study was to investigate MLN as the most promising candidate gene for LDA on proximal BTA23. We were able to detect a significantly LDA-associated single nucleotide polymorphism (SNP) affecting a predicted transcription factor binding site of the MLN gene which also changes the expression level of this gene. Therefore, this polymorphism is useful for a genetic test to reduce the incidence of LDA in GH cows.

Mutation analysis
First, comparative sequencing of the complete MLN gene including all exons and introns as well as the 39 and 59 UTR sequences was done in four GH cows affected by LDA and four control cows of the same breed. The genomic sequence of the control cows (FN298674) was not in complete agreement with the corresponding sequence of the Bos taurus genome 4.0 (NC_007324). In comparison to the reference sequence, the first intron of MLN was shorter by 557 bp in all GH cows sequenced. Eight German Fleckvieh cattle were also analysed for this sequence variation of MLN with the result of an identical length of intron 1 as the GH cows. The putative promoter of bovine MLN extends from 73 bp to 123 bp using the genomic MLN sequence identified here (FN298674). We detected a total of 32 polymorphisms, of which 30 were SNPs and two were short tandem repeats. All polymorphisms were tested for their information content in a sample of 48 GH cows, with one half of them affected by LDA. Only eleven polymorphisms were informative ( Table 1). None of the identified SNPs is located within the coding sequence of MLN, but FN298674:g.62G.A is located within the promoter and FN298674:g.90T.C within the non-coding first exon of MLN and the latter one affects a predicted NKX2-5 transcription factor binding site (Figure 1).
In order to correct for the data structure, we employed a logit model including a random sire effect besides the fixed genotypic SNP effects. This association analysis showed the most significant 2log 10 P-values for FN298674:g.90T.C and FN298674:g. 1891insG (3.8 and 3.7). The 2log 10 P-values of all other polymorphisms were distinctly smaller. In addition, we performed association analyses using a random additive genetic animal effect through a numerator relationship matrix and a heritability of h 2 = 0.3 for LDA estimated in the same population [1,2]. The most significant associations with LDA were found for FN298674:g.90T.C and FN298674:g.1891insG. The FN298 674:g.90T.C SNP showed the highest association with LDA. The SNPs FN298674:g.6689C.T and FN298674:g.6728G.A were no longer significantly associated with LDA. Therefore, we considered FN298674:g.90T.C and FN298674:g.1891insG as the most influential SNPs of MLN for developing LDA. The phenotypic variance explained after correction for the sire effect was at 3.1% for FN298674:g.90T.C and at 3.9% for both SNPs, FN298674:g.90T.C and FN298674:g.1891insG.

Haplotype analysis
The most common haplotypes for FN298674:g.90T.C and FN298674:g.1891insG comprising 95% of all cows were C-INS and T-WT. The haplotype including both mutant alleles (C-INS) was present in 55.1% of the LDA-affected and in 39.0% of the unaffected cows (Table S1). The haplotype trait analysis using these SNPs gave a x 2 value of 119.4 (p,0.0001) and x 2 values of 14.0 (p = 0.0002) to 60.8 (p,0.0001) for the distribution of the four single haplotypes among LDA-affected and control cows. Inclusion of further MLN polymorphisms did not improve the haplotype analysis.

Linkage analysis
We employed a multipoint non-parametric linkage analysis to re-assess linkage for MLN-associated markers in the data set previously used for the genome-wide search for LDA-QTL [4]. The chromosome-wide linkage analysis for 30 markers on BTA23 and all 14 paternal half-sib families showed a peak with a Zmean of 2.55 (P-value = 0.005) at the SNP FN298674:g.90T.C within MLN (Table S2, Figure S2).

Distribution of LDA-associated alleles in German Fleckvieh
The most significantly LDA-associated polymorphisms, FN298674:g.90T.C and FN298674:g.1891insG, were tested in 148 German Fleckvieh cattle. German Fleckvieh was used as a reference breed as this breed is known for its very low incidence of

Expression study
A qRT-PCR of cDNA gave evidence of a MLN expression in abomasal mucosa tissue. To test whether the FN298674:g.90T.C mutation might influence the MLN expression, samples of abomasal mucosa tissue from 55 previously genotyped GH cows were taken and analysed using qRT-PCR. Of these individuals, 26 were homozygous for the wildtype allele, 20 were heterozygous, and nine were homozygous for the mutant allele. The mean value of the DCT of the homozygous wildtype cows was used as calibrator [11]. The relative expression levels of MLN were decreased by 89% in cows being heterozygous or homozygous for the mutant allele C of the polymorphism FN298674:g.90T.C in comparison to cows homozygous for the wildtype allele T (Table 4). Among cows carrying the mutant allele C homozygously or heterozygously, the expression levels were identical. The error probability for the wildtype genotype versus the heterozygous or homozygous mutated genotype was at a P-value of 0.03 with a corresponding F-value of 4.8 (Table 4).

Discussion
This study demonstrated that the non-coding MLN transcription factor binding site mutation FN298674:g.90T.C is significantly associated with MLN expression levels and the incidence of LDA in GH cows. Furthermore, this polymorphism showed a significant linkage with LDA in GH cows from the same population. The MLN gene encodes the 22 amino-acid peptide motilin, a hormone which plays a role in stimulating gastrointestinal motility. Motilin is considered as an endocrine regulator and initiator of phase III of the migrating motor complex, the interdigestive peristaltic reflex in the gastrointestinal tract [5,12]. This hormone therefore causes gastric emptying in periods of fasting and its release ceases with the ingestion of food. Plasma levels of motilin increase cyclically every 90-120 minutes within the interdigestive periods, which evoke short time intervals of five  to ten minutes of strong peristaltic contractions initiating from the stomach [5]. In a study using the house musk screw as animal model for motilin studies, the contractile responses induced by motilin were dose-dependent [13]. There are also therapeutical approaches in using the positive effect of motilin agonists on gastric contractions concerning the treatment of diabetic gastroparesis [6,7] and altered motilin plasma levels are seen in some human diseases [14]. This study shows that MLN plays a role in the development of bovine LDA, a disease known to be initiated among other things by a reduced peristalsis of the gastrointestinal tract [3]. The reduced food-intake often seen in dairy cows after calving necessitates a high interdigestive peristalsis, which is a function of MLN. Naturally, its influence on the development of LDA is a fractional amount, because the disease is known to be genetically complex. After correcting the data for stratification, the polymorphism FN298674:g.90T.C within bovine MLN showed the highest association with LDA. This SNP affects a NKX2-5 transcription factor binding site. The effect of NKX2-5 on gene expression has been previously described [15] and SNPs within NKX2-5 binding sites have been reported to be functional for complex genetic diseases as systemic lupus erythematosus, rheumatoid arthritis, and graves disease in man as the mutations lead to a reduced expression of the respective gene [16]. Therefore, FN298674:g.90T.C is presumably the reason for the reduced expression of MLN in GH cows carrying a C allele at this position. Interestingly, while the association study showed the highest associations for the homozygous mutant genotype with LDA, the expression of MLN was reduced in cows with at least one copy of the mutated allele. This indicates a dominant effect of FN298674:g.90T.C. These results might be explained by the motilin-level regulation process. Motilin controls the gastrointestinal peristalsis reflex to provide gastric emptying during the interdigestive phases [5,12]. Plasma levels of motilin increase cyclically within these phases and evoke short time intervals of strong peristaltic contractions [5]. Therefore, though cows heterozygous for the FN298674:g.90T.C SNP show a considerably reduced basis expression of MLN, they might still be able to raise their motilin level more frequently or higher during the digestive periods where it is required and therefore reach a threshold to prevent LDA. Individuals homozygous for the mutated allele cannot raise their motilin level much beyond basis expression. Therefore, the gastric motility in heterozygous cows might be sufficient to prevent LDA.  The mean relative expression levels of MLN with their standard errors are given for the wildtype genotype T/T, the heterozygously (C/T) and the homozygously mutated genotype (C/C) standardized on the mean of the genotype T/T. Furthermore, the F-values and P-values for the differences in mean relative expression levels among the genotype T/T versus the genotypes C/T and C/C as well C/T versus C/C are given. RPL4 was used as housekeeping gene. doi:10.1371/journal.pone.0035562.t004 Though the homozygously mutated genotype of FN298674:g.90T.C significantly increases the risk of a cow to be affected by LDA and is the most likely causal SNP due to its effect on a transcription factor binding site, FN298674:g.1891insG within the first intron of MLN also showed a highly significant association with LDA. As FN298674:g.1891insG is in strong linkage disequilibrium (r 2 = 0.8) with FN298674:g.90T.C, its high association with LDA might be explained by this fact. However, the insertion mutation FN298674:g.1891insG might also contribute to the negative effect on MLN expression since an effect on a further regulatory domain can not be excluded.
To reveal the genotypic distribution of the SNPs FN298674:g.90T.C and FN298674:g.1891insG in a cattle breed with a low incidence of LDA [17], these SNPs were tested in German Fleckvieh. In German Fleckvieh, the wildtype allele was prevalent with a frequency of 0.653 at FN298674:g.90T.C and of 0.738 at FN298674:g.1891insG and the homozygous mutant genotypes for these polymorphisms showed an even lower frequency (0.161 and 0.114, respectively) than in GH cows unaffected by LDA (0.176 and 0.170, respectively). This corroborates the impact of MLN as a causal gene for LDA.
GH are high-performance dairy cows and due to their energy adapted-feeding and the cattle specific anatomic peculiarities they might be more susceptible for diseases caused by gastric motility decrease than other species. However, alongside cattle and human [18], gastric dilatation and displacement have also been reported in other species like dogs [19], cats [20], pigs [21], guinea pigs [22], and horses [23]. Therefore, this study might be of interest in the research of gastric motility disorders in these species.
In conclusion, we identified a SNP (FN298674:g.90T.C) which affects a predicted NKX2-5 binding site within the first noncoding exon of the bovine MLN gene. This polymorphism was shown to be correlated with a significantly lower expression of MLN and a significant increase of LDA incidence in GH cows. This is the first polymorphism showing a functional association with LDA and it can be used to test for LDA-susceptibility in this breed. MLN might be a useful candidate for the genetic analysis of gastric motility disorders in further species.

Ethics statement
All animal work has been conducted according to the national and international guidelines for animal welfare. All bloodsampling of LDA-affected cows was done in university clinics for cattle along with diagnostic protocols prior to surgery. The blood samples of unaffected cows were taken by veterinarians on farms. The project has been notified to the Lower Saxony state veterinary office, Niedersä chsisches Landesamt für Verbraucherschutz und Lebensmittelsicherheit, Oldenburg, Germany, and registered as a notifiable experiment with the number 33.9-42502-05-08A541.

Animals
We used EDTA-blood samples of a total of 1,136 GH cows. Of these animals, 601 cows had undergone a surgery because of LDA at the clinics for cattle of the veterinary universities Hannover or Giessen, Germany. The remaining 535 GH cows were unaffected by LDA. This sample was balanced by paternal grandsires. Each grandsire had to have at least one descendant affected by LDA and one unaffected descendant. The data set included 72 grandsires. Furthermore, EDTA-blood samples of 148 German Fleckvieh individuals with diverse ancestry were included into this study. Genomic DNA from EDTA-blood was extracted using the QIAamp 96 Spin Blood Kit (Qiagen, Hilden, Germany). In addition, samples from the abomasal mucosa tissue were taken from 55 GH cows at the slaughterhouse (Schlachthof Hannover and Schlachthof Giessen, Germany). All these animals were fasting prior to sampling. The samples were taken under equal conditions directly after slaughtering.

Sequencing and genotyping
The complete MLN gene was sequenced and scanned for polymorphisms in each four LDA-affected and unaffected GH cows. To ensure that the four control cows would not develop a LDA later in life, only cows with an age of at least nine years and without ever showing signs of LDA in their lifes were chosen for this first scan for polymorphisms. Sequence data were analysed using the Sequencher 4.9 software (GeneCodes, Ann Arbor, MI). PCR primers used for sequencing of MLN and for genotyping are given in Table 1 and Table S3. Sequencing of PCR products was performed on a MegaBACE 1,000 capillary sequencer (GE Healthcare, Freiburg, Germany). After this, each polymorphism identified was tested in 24 LDA-affected cows and 24 controls and, if it was informative, genotyped in a total of 1,136 GH cows. Markers were regarded as informative, if their PIC and heterozygosity were at 0.25 or higher. Genotyping of SNPs was performed using enzymatic digestion or sequencing (Table S4). For the genotyping of microsatellites, insertions, and deletions, each one primer was IRD-labeled and the PCR products were size-fractionated by gel electrophoresis on automated sequencers (LI-COR 4300, Lincoln, NE, USA) using 6% polyacrylamide denaturing gels (Table S4).

Expression analysis
For the analysis of MLN expression, RNA was isolated from abomasal mucosa tissue samples of 55 different cows using the NucleoSpinHRNA II kit (Macherey-Nagel, Düren, Germany) according to the manufacturer's instructions. The reverse transcription into cDNA was performed using 200 U SuperScript III Reverse Transcriptase (Invitrogen, Karlsruhe, Germany). A single PCR assay was used for quantification of the MLN gene transcript (GenBank NM_173938.2) using a forward primer situated on the boundary of exon 2 and exon 3 (59-TGC AGG AAA AGG AGA GGT ACA AG -39) and a reverse primer within exon 4 (59-GTT CAT CCT CAT TCC AAT TTC CA -39) amplifying an 155-bp product with a FAM-labeled TaqMan minor groove binding (MGB) probe (Applied Biosystems, Darmstadt, Germany) located at the boundary of exon 3 and 4 (59-CAG GAA GTT ATC AAG CTG A-39). Bovine RPL4 was used as a housekeeping gene. We used a forward primer within exon 2 (59-TTG CGC AAA AAC AAC AGA CA-39) and a reverse primer in exon 3 (59-GCC GGT ACC CCA AGA CTC A-39) amplifying a 81-bp product in combination with a VIClabeled TaqMan MGB probe (Applied Biosystems) spanning the boundary of exon 2 and 3 (59-TAG CAG GTC ATC AAA CC-39) according to the bovine reference mRNA sequence (GenBank NM_001014894). The quantitative real-time (qRT)-PCR was carried out using an ABI7300 sequence detection system (Applied Biosystems). The MLN transcript specific expression was normalised by the bovine RPL4 expression level (DCT), and the relative expression level was calculated by the DDCT method using the average DCT of the homozygous wildtype samples as calibrator [11]. All assays were performed in duplicates.

Statistical analyses
To test for association with LDA, case-control analyses based on x 2 -tests for genotypes and alleles were performed using SAS/ Genetics, version 9.3 (Statistical Analysis System, Cary, NC, USA, 2011). The ALLELE and HAPLOTYPE procedures of SAS/ Genetics were used for estimation of allele and haplotype frequencies, Hardy-Weinberg equilibrium and linkage disequilibrium among marker alleles. Logistic regression analysis was performed using the GLIMMIX procedure of SAS. The GLIMMIX procedure fits statistical models to data with correlations or nonconstant variability and where the response is not necessarily normally distributed and was therefore also used with our data. An animal model analysis was parameterized using a numerator relationship matrix for each individual and a heritability of 0.3 for LDA. This heritability was chosen as a mean value from former studies for LDA in GH cows where heritabilities of 0.2 to 0.5 were estimated [1,2]. Odds ratios (OR) and P-values were obtained from the CASECONTROL and GLIMMIX procedure. Promoter prediction was carried out by Promoter Prediction software (http://www.fruitfly.org). To detect the impact of intronic polymorphisms and to search for transcription factor binding sites, the TFSEARCH software (http://www.cbrc.jp/research/db/TFSEARCH.html) was used. Statistical calculation of pairwise linkage disequilibrium (LD) was performed using HAPLOVIEW 4.0 (http://www.broad.mit.edu/ mpg/haploview/). Statistical evaluation of MLN expression results was performed using standard procedures for comparing means with GLM (general linear model) of SAS, version 9.3.

Linkage Analysis
In addition to the association analysis, a multipoint nonparametric linkage analysis was conducted. For this purpose, 14 previously described paternal half-sib families segregating for LDA and containing a total of 360 animals were employed [4]. These families can be grouped into five grandsire families. We have taken the same microsatellites on BTA23 as in the previous study. In addition, we included the SNPs FN298674:g.90T.C and FN298674:g.1891insG within MLN as well as five new microsatellite markers located on the proximal end of BTA23. Positions of all markers were determined using a BLAST (http://blast.ncbi. nlm.nih.gov/Blast.cgi) analysis against the UMD_3.1 bovine genome assembly. The markers DIK5319 and BM47 could not be placed or were located on a different chromosome and therefore excluded from this study. A total of 30 markers on BTA23 were used. PCR was conducted according to a standard protocol. One primer of each pair was endlabeled with fluorescent IRD700 or IRD800. Primers and positions of the five newly developed microsatellites are given in Table S5. For the analysis of the microsatellite genotypes, the PCR products were sizefractionated by gel electrophoresis on an automated sequencer (LI-COR 4300, Lincoln, NE, USA) using 6% polyacrylamide denaturing gels. Multipoint non-parametric linkage analyses based on haplotypes identical-by-descent (IBD) for a half-sib design were employed [24,25]. The statistical analyses were performed using MERLIN, version 1.1.2 [26].