Associated regions for multiple birth in Brown Swiss and Original Braunvieh cattle on chromosomes 15 and 11.

Twin and multiple births have negative effects on the performance and health of cows and calves. To decipher the genetic architecture of this trait in the two Swiss Brown Swiss cattle populations, we performed various association analyses based on de-regressed breeding values. Genome-wide association analyses were executed using ~600 K imputed SNPs for the maternal multiple birth trait in ~3500 Original Braunvieh and ~7800 Brown Swiss animals. Significantly associated QTL were observed on different chromosomes for both breeds. We have identified on chromosome 11 a QTL that explains ~6% of the total genetic variance of the maternal multiple birth trait in Original Braunvieh. For the Brown Swiss breed, we have discovered a QTL on chromosome 15 that accounts for ~4% of the total genetic variance. For Original Braunvieh, subsequent haplotype analysis revealed a 90-kb window on chromosome 11 at 88 Mb, where a likely regulatory region is located close to the ID2 gene. In Brown Swiss, a 130-kb window at 75 Mb on chromosome 15 was identified. Analysis of whole-genome sequence data using linkage-disequilibrium estimation revealed possible causal variants for the identified QTL. A presumably regulatory variant in the non-coding 5' region of the ID2 gene was strongly associated with the haplotype for Original Braunvieh. In Brown Swiss, an intron variant in PRDM11, one 3' UTR variant in SYT13 and three intergenic variants 5' upstream of SYT13 were identified as candidate variants for the trait multiple birth maternal. In this study, we report for the first time QTL for the trait of multiple births in Original Braunvieh and Brown Swiss cattle. Moreover, our findings are another step towards a better understanding of the complex genetic architecture of this polygenic trait.

INTRODUCTION

The length of calving intervals and the number of calves born alive have a strong influence on the productivity of dairy cattle. Intensive selection, especially for milk yield in dairy cattle over the last 100 years, has led to a decline in female fertility due to negative genetic correlations between milk yield and female fertility (Miglior et al., 2017). To reverse the negative drifts in dairy cow fertility, producers benefited from the development of genomic selection and advanced management strategies; however, further efforts are needed.

Cattle are usually monoecious, so that a pregnancy usually ends with the birth of a single calf. Multiple births in cattle are rare. The multiple birth rate (MBR) varies between 1.02 and 9.6% depending on breed and study (Johanson et al., 2001; Moioli et al., 2017; Weller et al., 2008). They are generally lower in beef cattle than in dairy cattle (Atteneder, 2007; Lett & Kirkpatrick, 2018). Most of multiple births are due to multiple ovulations, as several ovulatory follicles mature at the same time. Around 5–10% of bovine twins are identical twins (Atteneder, 2007; Silva del Río et al., 2006). So far, however, the trait of multiple births has neither been studied in general nor with a genetic model for Brown Swiss (BS) and Original Braunvieh (OB) breeds.

Generally, multiple births are undesirable in dairy cows for several reasons. In particular, twin and multiple births are linked to increased health problems for the dam and calves. For example, there is a higher risk of ketosis, metabolic disorders, retained placenta, and displaced abomasum (Echternkamp & Gregory, 1999; Gregory et al., 1990; Pardon et al., 2012). In addition, the impact of multiple births on subsequent fertility was observed through their negative influence on conception rate and calving interval (Andreu‐Vázquez et al., 2012; Echternkamp & Gregory, 1999). Calves born in multiples have an increased risk for deficiency syndrome and mortality (Atteneder, 2007; Silva Del Río et al., 2007). Furthermore, a higher incidence of dystocia, abortion, and stillbirth has been found in multiple births (Atteneder, 2007; Gregory et al., 1990). In summary, these factors lead to higher production costs for the farmers. Consequently, selection against multiple births could improve fertility and profitability in dairy production.

Several non‐genetic factors have been described that could have an influence on MBR such as cow parity, and further environmental factors as the season of birth and the herd (Johanson et al., 2001; Silva Del Río et al., 2007). A significantly higher MBR was observed in multiparous cows (5.22–7.35%) compared to cows in first parity (1.63%) (Johanson et al., 2001). So far, two studies have analysed the relationship between milk yield and MBR. These studies presented contradictory results and thereby do not give a clear picture of the relationship between these two traits (Masuda et al., 2015; Murillo‐Barrantes et al., 2010). Most studies reported a higher MBR for births in summer months (Johanson et al., 2001; Silva Del Río et al., 2007). This suggests that the number of multiple ovulations is highest in late summer and autumn. Several studies indicate a positive phenotypic trend for MBR over time (Fitzgerald et al., 2014; Moioli et al., 2017; Silva Del Río et al., 2007), suggesting that MBR is correlated with further selection traits. Regarding the genetic trend, Ghavi Hossein‐Zadeh et al. (2009) found a negative trend and Murillo‐Barrantes et al. (2010) found no change in the genetic component of MBR over time. In most dairy cattle populations, MBR seems to be increasing at the phenotypic level. The heritability estimates for multiple births range from 0.011 to 0.160, suggesting that a small but non‐zero proportion of the observed trait variation can be attributed to genetic factors (Fitzgerald et al., 2014; Johanson et al., 2001; Lett & Kirkpatrick, 2018). Estimated values using a linear model were lower than those using threshold models. Furthermore, an Austrian study reported higher estimates for the heritability of the maternal multiple birth trait (0.017–0.063) compared to the estimates for the heritability of the direct multiple birth trait (0.001–0.005) (Atteneder, 2007).

Maternal multiple birth has been studied with different QTL mapping strategies in Norwegian cattle, as well as in Israelian or North American Holsteins using family‐based microsatellite interval mapping approaches (Bierman et al., 2010; Lien et al., 2000; Weller et al., 2008). They detected multiple QTL on 13 different chromosomes, depending on the study conducted and the population examined (Cobanoglu et al., 2005; Lien et al., 2000; Weller et al., 2008). Paternal half‐sib families for North American Holstein cattle were analysed using single‐marker models or combined linkage–linkage disequilibrium approaches, which resulted in the discovery of multiple QTL on eight chromosomes (Bierman et al., 2010; Kim et al., 2009). The Italian Maremmana beef breed was analysed in a study using a single‐trait linear mixed effect model as well as an animal threshold model including the number of calves born per cow as the phenotype. They identified a single significantly associated SNP on chromosome (chr) 24 by a genome‐wide association studies (GWAS) based on 54‐k SNP data (Moioli et al., 2017). So far, a study has analysed the Holstein population using large‐scale data and identified a major QTL for multiple births on chr 11 with a possible effect on the genes LHCGR and FSHR (Widmer et al., 2021). This QTL region was confirmed in a recent study in the North American Holstein population using whole‐genome sequencing (WGS) data (Lett & Kirkpatrick, 2022). It is therefore very likely that multiple births are a polygenic trait where almost no QTL have been identified yet. The study presented in this paper is based on a large‐scale dataset of phenotypes and genotypes for the current Swiss dairy cattle breeding populations of BS and OB.

We aimed to perform an extensive genetic analysis of the multiple birth trait in Swiss OB and BS cattle. To predict breeding values, phenotypic data for single and multiple births from 2 decades were used, which were recorded in the national animal database, as well as data from pedigrees. In addition, extensive genotype data obtained from routine genomic selection of male and female animals were used to detect associated QTL using GWAS and haplotype regression analysis. To identify linked genomic regions and candidate causal variants, different fine‐mapping procedures were performed. In a final step, we evaluated the effect of the identified haplotype on the routinely recorded birth and fertility traits.

MATERIALS AND METHODS

Phenotypes

Phenotypic data were available on a large scale through the Swiss national animal recording database. In this study we used data collected between 2006 and 2018 from the breeding organisation Braunvieh Schweiz (Zug, Switzerland). The raw dataset contained 5 692 022 birth records including the multiple birth code from BS and OB cattle. Data preparation and analysis for variance component estimation and breeding value estimation were done with internal software written in R (R Core Team, 2020). Birth records subsequent to an embryo transfer were excluded. For the genetic studies of the discrete multiple birth trait 1 367 529 records remained after data preparation and validation. The overall MBR was 7.83%. Further information regarding the final dataset is shown in Table 1.

TABLE 1

Estimated effects of the factor levels for the fixed effects of the estimated breeding values based on the final dataset

Fixed factor	Level	Number of observations per level a	Estimated effect of factor level
Parity	1	351 678	0.212
	2	308 988	0.271
	3	222 356	0.286
	4	168 837	0.295
	5+	315 670	0.302
Use of sexed semen	No	1 312 571	0.269
	Yes	54 958	0.264
Season of birth	Spring	231 908	0.261
	Summer	266 155	0.278
	Fall	519 258	0.270
	Winter	350 208	0.265
Number of inseminations	1	876 021	0.270
	2	324 011	0.267
	3	106 625	0.266
	4	36 382	0.265
	5	13 697	0.256
	6	57 776	0.251
	7+	5017	0.253

Based on the final dataset of 1 367 529 records.

Variance components estimation

We fitted the following mixed linear model to our phenotypic data described above according to a previous study in a different breed (Widmer et al., 2021): where b is the vector of the fixed effects, h represents the vector of the random herd‐year effect, mbd and mbm are the direct (calf) and the maternal (dam) genetic effects respectively, and ε is the residual. The fixed effects of season, parity, use of sexed semen, and number of inseminations leading to pregnancy were taken into account. In Table 1, the numbers of observations per level of each fixed effect are illustrated. X, W, Z d, and Z m represent the design matrices for the random (W, Z d, and Z m) and fixed (X) effects. The selection of the fixed effects is based on previous unpublished work that analysed a similar dataset in Switzerland. The items in the vector of observations (y) were coded with 1 for single and 2 for multiple birth (twin or triplet). Dataset filtering was performed to eliminate all records from herds with fewer than 260 records per herd in total (<20 per year on average), and all records in herd‐year classes with fewer than five records. Finally, this resulted in a dataset of 419 111 records. Estimation of variance components were conducted for BS and OB together by using the software vce (Neumaier & Groeneveld, 1998).

Breeding value prediction

Breeding values prediction was performed with the MiX99 software (MiX99 Development Team, 2017; https://www.luke.fi/mix99) by applying the mixed linear effects model shown above and run for BS and OB together. Within this analysis, the dataset filtering excluded records from herds with fewer than 260 records but without defining a minimal number of records per herd‐year levels. In total, 1 367 529 records were available for breeding value prediction. The estimates from the prior step were the base for the required values of variances and covariances. The estimation of the reliability of all breeding values was based on an approach by Tier and Meyer (2004) using the apax99 program in the MiX99 software. Furthermore, breeding values were standardised to a mean of 100 and a standard deviation of 12.

The de‐regression of the estimated breeding values (EBVs) for the direct (mbd) and maternal (mbm) multiple birth traits was performed according to Garrick et al. (2009) (Table 2). For the association analyses, all animals having a de‐regressed EBV with a corresponding reliability above 0.35 were selected. The EBVs were not weighted by the reliability for the further steps. The following analyses were carried out separately for the BS and the OB breeds. They were separated using the pedigree‐based gene proportion. Minimum required gene proportion was defined at 0.7 for each of the tail populations. Finally, for BS 3792 and 7847 animals remained in the analysis for the traits mbd and mbm respectively. For OB there were 1332 and 3508 animals used for the following analyses for the mbd and mbm traits.

TABLE 2

Summary statistics of the de‐regressed breeding values for the direct and maternal multiple birth traits

Trait	Min	Max	Mean	SD	n observations
mbd	−50.810	72.971	−0.065	6.441	3293
mbm	−55.453	186.960	0.892	9.268	7284

Abbreviations: mbd, multiple birth direct; mbm, multiple birth maternal.

Genotypes

For 97‐k animals, routine SNP genotype were available. Imputation was performed as described in a previous study (Widmer et al., 2021) for the SNP data generated for genomic selection. The reference dataset for the 150‐k array included 6370 BS and 1516 OB animals and for the high‐density array 1686 BS and 421 OB cattle. The final marker set contained 110 510 and 681 178 SNPs for each density (150 k and high density) respectively. SNPs were filtered separately for the BS and OB populations using the following thresholds: minor allele frequency >0.01 and SNP call rate >0.99 in the genotypic data of the reference population. The ASR‐UCD1.2 bovine assembly was used as the reference genome for the SNP data.

Association studies

Single SNP regression

Genome‐wide single SNP association studies were performed using a mixed model in the software snp1101 (Sargolzaei, 2021). The genomic relationship for the cattle used in the analyses was calculated in order to correct for population stratification in the model (VanRaden, 2008). The aim of this approach was to detect variants that were significantly associated with the studied traits mbd and mbm.

Bayes B approach

To fit one genotype after another can simply lead to biases due to LD and stratification (Fernando et al., 2017). Therefore, fitting subsets of genotypes at the same time can address this paradigm. Window‐based association analyses were carried out using GenSel software (Putz, 2021) and the BayesB algorithm (Fernando & Garrick, 2013). First, the proportion of loci with zero effect (known as π parameter) was estimated from the same dataset, and a value of 0.989 was set as the a priori starting value. We estimate the window variance for genomic windows of 1‐Mb length.

Haplotype analysis

We used a haplotype regression analysis restricted to the previously identified genomic windows, explaining a significant proportion of the total genetic variance, to detect the haplotype which include the presumed causal QTL. Following a previously described approach by Pausch et al. (2016), we aimed to identify significantly associated haplotypes for both breeds within each window. We intersected haplotypes using all odd numbers between 9 and 301 SNPs as lengths within the identified segment per breed. We slid the starting point of the haplotype SNP‐wise and used different lengths so that we can consider each haplotype for the effect estimation. The de‐regressed EBV were again used as response variables for the estimation of haplotype effects. For the following fine‐mapping approach to identify any candidate causal variants using whole‐genome sequence data, we selected the most significantly associated haplotype (lowest p‐value).

The top‐associated haplotype was analysed regarding associations with routinely recorded and available calving and fertility traits. Therefore, the haplotype was fitted as a single diplotype representing haplotype carrier states (0 = non‐carrier, 1 = carrier and 2 = homozygous carrier). To correct for population stratification, the genomic relationship was included in the model in the GCTA software package (Yang et al., 2011). Additionally, we performed an analysis to investigate a possible co‐association between the most significantly associated haplotype for our trait of interest in each breed and further recorded fertility traits. For this test we used routinely available de‐regressed EBV for all fertility and birth traits.

Fine‐mapping

WGS data were searched for variants in significant LD with the selected most significantly associated haplotype, which was decoded as diplotype using the plink software package (Purcell et al., 2007). For BS, the 1‐Mb window between 75 and 76 Mb on chr 15 was selected based on the associated Bayes B window and the localisation of the highly associated haplotypes. For the OB breed, we selected the 1‐Mb window between 88 and 89 Mb on chr 11. WGS data provided by the 1000 Bull Genomes Project run 8 were used. This dataset contains 4109 animals and includes 237 BS and 81 OB cattle (Hayes & Daetwyler, 2019). We converted the data chromosome‐wise from VCF format to plink input files using the software VCFtools v0.1.16 (Danecek et al., 2011). Subsequently, we transformed the data into binary files and performed VCF quality controls using the plink v.1.90b3.46 software (Purcell et al., 2007). We removed variants with a missing call rate >0.1, a minor allele frequency <0.01, and samples with missing call rate >0.1. In addition. Samples and variants with a Mendelian error rate >2% were excluded. For the animals of the 1000 Bull Genomes Project, pedigree data were provided by the Braunvieh Schweiz (Zug, Switzerland) and included 526 duos and 19 trios that could be used for Mendelian analyses. Of the 237 whole‐genome sequenced BS and 81 OB animals, 137 samples were used for LD‐based variant detection after quality control. Only for these animals could additional SNP genotype data be obtained: 5933 variants were located in the 1‐Mb window on chr 15 for BS and, for the OB breed, the final analysis included 8041 variants in the 1‐Mb region on chr 11.

Furthermore, we have looked at the frequency of the top‐linked variants across breeds. We used all 4109 samples in the data file provided by the 1000 Bull Genomes Project. We analysed all breeds that contained at least 50 records.

RESULTS

Prediction of breeding values

The raw values resulting from the estimation of the variance components are shown in Table 3. The heritability of the direct (animal) and the maternal (dam) genetic effects was 0.2 10−4 and 0.040 respectively.

TABLE 3

Estimated raw variance components

Component	Raw value	Standard error
Herd‐year	0.88 × 10−4	0.15 × 10−4
Animal	0.72 × 10−6	0.18 × 10−5
Dam	0.16 × 10−2	0.74 × 10−4
Residual	0.037	0.91 × 10−4
Correlation animal/dam	−0.34 × 10−4	0.43 × 10−4

We estimated breeding values for the traits mbm and mbd. These values were used in a de‐regressed alteration as input variable for the association analysis. The estimated effects of the factor levels for the fixed effect of use of sexed semen, parity, season of birth, and number of inseminations leading to pregnancy are shown in Table 1. The use of sexed semen has a significant negative effect on the occurrence of multiple births. In addition, the prevalence of multiple births was lower in primiparous cows and in the summer months, the highest incidence for multiple births was observed. The empirical distributions of the raw, de‐regressed, and standardised values of the breeding value prediction can be seen as boxplots in Figure 1.

FIGURE 1

Raw, standardised, and de‐regressed breeding values for the traits direct (mbd) and maternal (mbm) multiple birth

The reliabilities of the standardised breeding values of 1 905 641 animals had a mean of 0.263 for mbd and 0.264 for mbm, and the standard deviation was 81.663 for mbd and 8.858 for mbm. As a measure of the genetic trend, we grouped the mean estimated breeding values by year of birth of the animals for each trait. Interestingly, no clearly discernible genetic trend was found for the trait mbm in either breed or for the trait mbd in BS (Figure 2). For the trait mbd in OB only, a negative trend was observed for 2000–2010.

FIGURE 2

Genetic trend of the estimated breeding values of multiple births from 2000 to 2019

A QTL close to the gene on chromosome 11 for OB

For OB cattle, the window‐based BayesB approach revealed a QTL with a significantly associated window that explained 5.88% of the genetic variance of the trait mbm (Figure 3a). The window is located on chr 11 between 88 000 459 and 88 996 149 bp (Table 4). The values of p > 0 (proportion of models where this window was included, and thus explained for >0% of the genetic variance) and of p > average (proportion of models where this window explains for more than the amount of variance that would be explained if each window had the same effect) were 1 and 0.980 respectively. Therefore, this QTL reached a significance threshold of 5% in the BayesB approach. The other three possibly associated windows, each located on different chromosomes, explained smaller amounts of genetic variance (<2.5%), and were identified below the significance level of 5% (Table 4).

FIGURE 3

Manhattan plot of genome‐wide association studies for the trait mbm. (a, c) results of the window‐based BayesB approach and (b, d) the single SNP regression with Bonferroni‐corrected threshold level of 5% (blue line) and 1% (red line). (a, b) Original Braunvieh, (c, d) Brown Swiss

TABLE 4

Associated genome regions from the BayesB window approach for the trait maternal multiple birth

Breed	Chr	Start position	End position	n SNP	Proportion of explained genetic variance (%)	p > 0	p > average
Original Braunvieh	11	88 000 459	88 996 149	366	5.88	1	0.980
	22	24 002 842	24 995 229	224	2.36	0.985	0.754
	19	36 007 788	36 998 151	274	1.51	0.967	0.729
	1	107 022 412	107 999 432	304	1.00	0.980	0.526
Brown Swiss	15	75 001 650	75 995 887	339	3.77	0.997	0.950
	5	57 008 008	57 996 196	196	2.13	0.962	0.664
	11	24 001 475	24 994 135	296	1.83	0.995	0.897

Abbreviations: Chr, chromosome; p > 0, proportion of models where this window was included, and thus explained for more than 0% of the genetic variance; p > average, proportion of models where this window explains for more than the amount of variance that would be explained if each window had the same effect.

In the single SNP regression models for the trait mbm, no significantly associated SNPs were found at the Bonferroni‐corrected threshold of 5% (Figure 3b). Different suggestive signals can be observed on chr 1, 10, 11 and 22. The significantly associated window on chr 11 from the BayesB analysis overlapped with the suggestive QTL on chr 11. The top SNP is located at position 88 765 206 bp. For the second trait analysed (mbd), we could not observe a significant signal in either of the two GWAS analyses.

We selected a 1‐Mb‐segment on chr 11 at 88 Mb for haplotype association analyses based on the results from the BayesB analysis shown above. The 16 top‐associated haplotypes are listed according to their p‐value from the association analysis in Table S1. They were located in the region between 88.74 and 88.83 Mb on chr 11 and contained 11 to 15 SNPs. We selected from the most significantly associated haplotypes the longest haplotype (chr11: 88 748 439–88 808 495) for further analysis. This haplotype comprises 15 SNPs and had a frequency of 0.583 in the OB genotype dataset. For this haplotype, we discovered an estimated additive effect of −3.573 (± 0.751 standard error), indicating a negative effect on the trait of interest mbm.

LD‐analysis revealed six variants located between 88.75 and 88.79 Mb on chr 11 showing r 2 values ≥0.7 with the chr11: 88 748 439–88 808 495 haplotype (Table S1), but none were in perfect LD (r 2 = 1; Figure 4b). By far the highest r 2 value (0.905) was determined for the variant chr11: 88 791 842 A>T (Figure 4b), which, like the other five variants, mapped to a non‐coding region on chr 11 (Table 5). The closest annotated gene is ID2, which codes for the inhibitor of DNA binding 2 and is located between 88 604 880 and 88 607 401 bp (Figure 4a). All these six variants are located within the top‐associated haplotype between 88.75 and 88.81 Mb as described above (Table S1).

FIGURE 4

Graphical representation of the associated QTL region on bovine chromosome 11 and linkage disequilibrium (LD) analysis results for the mbm trait in Original Braunvieh. (a) Screenshot of the region between 88 and 89 mb from www.ensembl.org including the localisation of the genes in the region. The blue bar shows the localisation of the top‐associated haplotype. The red star represents the variant with the highest LD score. (b) The variants in r 2 ≥ 0.7 with the top‐associated haplotype are highlighted and shown with their genomic position. (c) Heatmap showing the LD in the region between 88 and 89 mb highlighting the haplotype that was added as an additional variant

TABLE 5

Overview of the highly associated variants from the linkage disequilibrium analysis for the trait multiple birth maternal

Breed	Chr	Position a	Variant a	Impact	Associated gene a	MAF	LD (r 2)
Original Braunvieh	11	88 791 842	A>T	Intergenic variant	5′ of ID2	0.343	0.905
Brown Swiss	15	75 213 046	T>C	Intron variant	PRDM11	0.142	0.970
		75 297 912	C>T	3′ UTR variant	SYT13	0.142	0.970
		75 399 114	G>T	Intergenic variant	5′ of SYT13	0.142	0.970
		75 402 900	G>A	Intergenic variant	5′ of SYT13	0.140	0.970
		75 405 408	T>G	Intergenic variant	5′ of SYT13	0.142	0.970

Abbreviations: Chr, chromosome; UTR, untranslated region.

ASR‐UCD1.2/bosTau9 assembly.

Since candidate causative variants for the identical trait may differ between breeds of the same species, we have looked at the frequency of the highly linked variant on chr 11 across breeds. The frequency of this top‐linked variant (chr11: 88791842 A>T) was high in several of the breeds including BS (Figure S1a). It is obvious that this variant segregates strongly in breeds other than OB, even in unrelated breeds.

A QTL close to the and genes on chromosome 15 for BS

For the studied BS population, the window‐based BayesB approach for BS revealed evidence for a QTL with a significantly associated window that explained 3.77% of the genetic variance of the trait mbm (Figure 3c). The localisation of this window is on chr 15 between positions 75 001 650 and 75 995 887 bp (Table 4). The values p > 0 and p > average reached values of 0.997 and 0.950 respectively. Therefore, the QTL from the BayesB analysis was significant at a threshold of 5%. Two other associated windows, each located on different chromosomes, were identified below the significance level of 5% (Table 4). They explained <2.5% of the genetic variance.

The result from single SNP regression analysis did not show any significantly associated SNPs at the Bonferroni corrected level of 5% for the mbm trait (Figure 3d). Nevertheless, a suggestive signal was observed on chr 15. The best‐associated SNP is at position 75 847 291 bp. This overlapped with the significantly associated window on chr 15 from the BayesB analysis. Further suggestive but non‐significant regions can be identified on chr 2 and 6 from the single SNP GWAS models. Also, for the BS population, we could not find any significant association signals for the trait mbd in either of the two GWAS analyses.

Based on the results from the BayesB analysis shown above, we used a 1‐Mb segment on chr 15 starting at 75 Mb for haplotype association analyses. The six best‐associated haplotypes according to their significance level from the association test were located in the region between 75.43 and 75.56 Mb on chr 15 and contained 41–51 SNPs (Table S1). The longest top‐associated haplotype (chr15: 75430142–75 561 032) was selected for further analysis, for which we identified an estimated additive effect of 5.018 (±0.686 standard error), which indicates that it has a positive effect on the trait mbm. This haplotype included 51 SNPs and had a frequency of 0.583 in the entire genotype dataset of the BS population studied.

The analysis of the pairwise LD (r 2) between the top‐associated haplotype (chr15: 75 430 142–75 561 032) and all detected variants in the 1‐Mb region revealed 60 variants located between 75.09 and 75.96 Mb on chr 15 with r 2 values ≥0.7 (Table S1). Five out of these 60 variants had values of r 2 ≥ 0.95 (Table 5), but for none was perfect LD (r 2 = 1) observed with the haplotype (Figure 5b). Of these five highly associated variants, one maps in an intron of the PRDM11 gene, one is located in the 3′ UTR of the SYT13 gene, and three are intergenic variants located 5′ upstream of the SYT13 gene (Table 5 and Figure 5a). None of the highly associated variants (r 2 ≥ 0.95) and only 6 of these 36 variants with r 2 values ≥0.7 map into the region of the top‐associated haplotype between 75.43 and 75.56 Mb (Table 5 and Table S1).

FIGURE 5

Graphical representation of the associated QTL region on bovine chromosome 15 and linkage disequilibrium (LD) analysis results for the mbm trait Brown Swiss. (a) Screenshot of the region between 75 and 76 mb from www.ensembl.org including the localisation of the genes in the region. The blue bar shows the localisation of the top‐associated haplotype. The red stars represent the variants with the highest LD scores (r 2 ≥ 0.95). (b). The variants in r 2 ≥ 0.7 with the top‐associated haplotype are highlighted and shown with their genomic position. (c) Heatmap showing the LD in the region between 75 and 76 mb highlighting the haplotype that was added as an additional variant

Additionally, we analysed the occurrence of the five chr 15 variants that were in high LD (r 2 ≥ 0.95) with the top‐associated haplotype across breeds. Interestingly, the frequencies of these five variants were low (<0.1) in all breeds except BS (Figure S1b–f). Four out of five variants segregated only in a small number of breeds. Only the variant chr15: 75 405 408 T>G occurred in all analysed breeds (at low level in OB), with the only exception of Angus where the variant does not segregate.

No detectable haplotype effects on routinely available fertility and calving traits

We estimated effects of the haplotypes for each breed, OB and BS, on routinely accessible fertility and calving traits to analyse a potential co‐association between other traits and the most significantly associated haplotypes for mbm. For the top‐associated haplotypes in OB (chr11: 88748439–88 808 495) and BS (chr15: 75430142–75 561 032) we found no significant association at the Bonferroni corrected level of 5% (p‐value ≤0.038; Table S1).

DISCUSSION

Impaired female fertility is a well‐known problem in high‐performance dairy cattle and multiple births have been an undesired trait for a long time. In this study, based on large‐scale genotyping and phenotyping data, we identified genetic factors that are highly likely to affect this trait. Using the de‐regressed breeding values for the trait of direct and maternal multiple births as phenotypes, two novel QTL were detected in two breeds: a QTL on chr 11 in OB and a QTL on chr 15 in BS. For OB, we found a variant located in the 5′‐regulatory region of the ID2 gene as potential candidate causal variant. In the BS population, we detected five variants that could presumably affect the expression of the PRDM11 and SYT13 genes and therefore could potentially be causal candidates.

The genetic parameter estimates (variance components of direct and maternal genetic effect for trait multiple birth) were similar to those reported before in different cattle populations (Ghavi Hossein‐Zadeh et al., 2009; Weller et al., 2008; Widmer et al., 2021). Since the phenotypic observations do not have a normal distribution, it is non‐trivial to develop a method to predict breeding values for a binary trait such as multiple births. However, the use of linear mixed models can lead to sufficient rankings of animals according to predicted breeding values, as Negussie et al. (2008) have shown, but the results obtained may not be valid and predictions do not have the same evidence as the discrete response variable. These challenges could be the reason for the extreme outliers in the standardised breeding values that we observed for both traits analysed (mbd and mbm). The potential of a threshold model or a generalised linear mixed model in the context of this trait needs further research. Previous studies have already used threshold models or generalised linear mixed models for genetic evaluations and predictions in a general purpose and for the analysis of multiple births (Johanson et al., 2001; McGovern et al., 2021).

In the last 2 decades, no clear genetic trend could be observed for both studied traits (mbd and mbm) in BS and for the trait mbm in OB. OB data for the trait mbd showed a negative trend from 2000 to 2010, in contrast to earlier results in Norwegian cattle where a positive trend was observed (Karlsen et al., 2000). Since the reliability and heritability of the trait mbd are low, the genetic trend could be caused by non‐genetic effects. However, it is very likely that multiple births were generally not selected in the Swiss populations studied, neither by direct selection nor by genetic correlation with selected traits. Thus, the development of breeding value estimation could be an important selection tool to reduce the risk of multiple births, hopefully leading to improved animal health and welfare of dairy cattle. Implementing a selection program into the local breeding schemes could be the next step.

In our datasets, we found significant associations for the maternal trait. This shows the important role of the maternal component for the complex trait of multiple births. The QTL on chr 11 for OB and on chr 15 for BS were detected by the window‐based BayesB approach. Only suggestive signals were observed in the single SNP regression analysis. As shown in this study, merely a small part of the genetic variance of quantitative traits can be significantly revealed by the identification of a single marker (Maher, 2008; Visscher et al., 2010). The BayesB approach considers all markers simultaneously as random effects and can therefore explain most of the genetic variance (Fan et al., 2011; Hayes et al., 2010). A window‐based approach captures the majority of the variability at an associated trait locus (Fernando & Garrick, 2013).

The identified QTL for OB explained 5.88% and the QTL for BS 3.77% of the genetic variance of the trait. Because mbm is a classical polygenic trait, a larger dataset could identify additional QTL with smaller effects. The non‐significantly associated signals resulting from the window‐based BayesB analysis on three additional chromosomes for OB and two additional chromosomes for BS (Table 4) may be suggestive QTL that need to be investigated in future studies. The theory of a polygenic trait is reinforced by previous studies showing multiple QTL on different chromosomes (Bierman et al., 2010; Kim et al., 2009; Weller et al., 2008; Widmer et al., 2021). Regarding chr 11 and 15, only combined linkage–linkage disequilibrium analysis for North American Holstein sires and ANOVA analysis in Israeli Holsteins showed QTL for twinning rate on the same chromosomes but located in different segments (Kim et al., 2009; Weller et al., 2008). Interestingly, no QTL for the multiple birth trait is reported in the bovine QTL database in the genome regions described in this study (Cattle QTL Database, 2022). There is little evidence in the literature of a major QTL for the multiple birth trait in cattle. In an own study, a major QTL on chr 11 at 31 Mb was detected with a candidate causal variant in the 5′ regulatory region of the LHCGR gene in Swiss Holstein (Widmer et al., 2021) and recently confirmed in the North American Holstein population (Lett & Kirkpatrick, 2022). In two older studies with data from Norwegian cattle and North American Holsteins, a candidate gene with a possible major effect on the trait, IGF1, was observed (Kim et al., 2009; Lien et al., 2000). The region on chr 5 encompassing IGF1 does not intersect with the suggestive QTL region identified for BS using the BayesB approach in our analysis. For bovine ovulation rate, a single study provided strong evidence for a major QTL in the region of 13.6 to 14.8 Mb on chr 10 containing three possible candidate genes (SMAD3, SMAD6, and IQCH) (Kirkpatrick & Morris, 2015). Presumably, a global dataset could improve our analyses and uncover additional QTL for the trait of interest.

LD‐based filtering of sequence variants located in the QTL regions did not reveal a single variant that had perfect LD with the top‐associated haplotypes for either breed. However, for BS on chr 15, a total of 60 variants showed high scores of LD with the haplotype (r 2 ≥ 0.7) forming an 870‐kb block. This long range is to be expected due to the previously described relatively long‐range LD pattern in cattle populations (de Roos et al., 2008). As previously shown, a single genomic region can harbour multiple QTL for a polygenic trait (Thaller et al., 2003); however, it is more likely that the identified long segment is due to LD between adjoining variants and only one of them represents the actual causal variant. It could be suspected that one of the five variants in highest LD (r 2 ≥ 0.95) is the most likely one (Table 5). One of those variants is a variant in the intron of the PRDM11 gene located between 75 202 012 and 75 285 991 bp on the forward strand on chr 15. This gene codes for the PR‐domain containing protein 11 that is predicted to enable chromatin binding activity It is involved in several processes, including negative regulation of cell growth, positive regulation of the apoptotic process in fibroblasts and regulation of transcription, but no role in reproduction has yet been reported (Fog et al., 2015). Different members of the PRDM family have an influence on germ cell specification with PRDM9 being expressed during meiosis (Fog et al., 2012). The other four variants with r 2 ≥ 0.95 are close to the SYT13 gene, which is located on chr 15 between 75 297 582 and 75 339 698 bp on the reverse strand, and may have a regulatory influence on its expression. SYT13 codes for synaptotagmin 13, which belongs to the large family of synaptotagmin proteins whose members function as membrane transporters in multicellular organisms (Quiñones‐Frías & Littleton, 2021). SYT9, a member of this family, affects the release of follicle stimulating hormone (FSH) in female mice and thus the oestrus cycle and ovulation (Roper et al., 2015). Therefore, we postulate that these five variants are possible candidates for a causal effect on the trait multiple births in BS cattle.

For the OB population, six variants with higher LD (r 2 ≥ 0.7) were observed in a 40‐kb window on chr 11. By far the highest r 2 with 0.905 was observed with the variant chr11: 88791842 A>T, a variant close to the ID2 gene located between 88 604 880 and 88 607 401 bp on the reverse strand, which encodes for the inhibitor of DNA binding 2. This gene belongs to the inhibitor of DNA binding family whose members are transcriptional regulators. Despite the nearly 200 kb distance to this gene, this non‐coding variant could probably have a regulatory effect on ID2 expression. Such a possible impact on ID2 is supported by the block of increased LD in this region, as shown by the heatmap presented (Figure 4c). ID2 is expressed at different levels during the oestrus cycle in the ovarian tissue in mice and may play a role in negatively regulating cell differentiation (Zavareh et al., 2018). In cattle, the expression of maternal ID2 transcripts decreased from immature to mature oocyte (Thélie et al., 2007). Expression of ID2 in pigs is influenced by FSH and cumulus oocyte complexes (Verbraak et al., 2011). In birds, ID2 is sufficient for the expression of the follicle stimulating hormone receptor (FSHR) (Johnson & Woods, 2009). Moreover, the expression of ID2 is significantly decreasing in granulosa cells of pre‐ovulatory follicles compared to mid‐oestrous follicles in mares (da Silveira et al., 2014). Therefore, we postulate the possibility of an influence of a regulatory variant of ID2 on ovulation and thus on the occurrence of multiple births. The candidate variant (chr11: 88791842 A>T) segregates in various breeds at high frequency (Figure S1a). Therefore, it could be assumed to represent an old mutation that probably arose before the formation of modern breeds. A validation analysis was performed and we were not able to confirm an association of that specific variant in the Swiss BS, Simmental or Holstein cattle populations (data not shown). Therefore, it is more likely that another variant in the QTL region that is in high LD with chr11: 88791842 A>T is causally effective, or that multiple variants have the same effect on the trait.

Since the proposed candidate causal variants for multiple births in the BS breed (Figure S1b–f) segregate in various breeds at low frequency, they probably occurred before the formation of modern breeds. A validation analysis as described for OB was not possible, because for most of the breeds analysed, neither phenotypes nor genotypes were available to investigate a possible effect on the trait multiple births. Remarkably, none of the five variants that had an LD score >0.95 belonged to the interval of the top‐associated haplotype for BS, as they mapped further upstream on chr 15. Interestingly, the haplotype and the variants are not in a block of raised LD as shown in the presented heatmap (Figure 5c), which could explain this phenomenon. The reported variants found in this study strongly suggest a possible functional connection to the genes SYT13 and PRDM11 in BS and ID2 in OB cattle respectively.

The observed associated effect of the specified haplotype on chr 11 on the investigated trait mbm was negative in OB cattle. Female haplotype carriers can therefore be expected to have a lower incidence of multiple births. By contrast, the effect for the haplotype on chr 15 for BS was positive. Consequently, carriers have a higher risk for multiple births. The assessment of the potential impact of haplotypes on other available birth and calving traits showed no significant impact.

CONCLUSIONS

By analysing large‐scale genotype and phenotype data in the two Swiss Braunvieh cattle populations, we have identified independent QTL for maternal multiple birth on chr 11 for OB and on chr 15 for BS. Furthermore, we provide evidence for linked variants in both breeds that possibly affect expression of ID2 as well as STY13 and PRDM11 that might impact the occurrence of multiple births. Therefore, these findings improve the understanding of the genetic architecture of this complex trait in cattle and for mammalian reproduction in general. Further studies are desirable, in particular to provide functional evidence of the causal relationship we postulated between the QTL discovered and the phenotype studied.

CONFLICT OF INTEREST

None declared.

Figure S1

Table S1

Table S2

Table S3

Table S4

Table S5

Click here for additional data file.