Implications
Twin birth is disadvantageous for most beef and dairy producers, because of its association with several unfavourable effects. However, in economically marginal production systems on range conditions and seasonal calving, a local suckler breed with high resilience, excellent calving and good mothering abilities could be proficiently selected for increasing twinning rate. By weaning two calves, yearly meat outcome per cow would increase. This work aims to investigate the feasibility of selection for twinning in the local Italian breed of Maremmana cattle.
Introduction
Every year, Italy imports about one million cattle; in 2011, 1 341 000 animals were imported, mainly for fattening (FAO, 2014). Improving fertility through an increase of twinning rate in beef herds could help in relieving the scarcity of calves. However, negative consequences of twinning such as dystocia, calf abandoning or calf mortality need to be evaluated in each specific environment and breed. Twinning occurs at a low frequency in cattle; according to Vinet et al. (Reference Vinet, Drouilhet, Bodin, Mulsant, Fabre and Phocas2012), twinning rate (number of twin births/100 births) ranges from 0.6% (in American Hereford) to 5% (in French Maine-Anjou and Holstein). Twinning in cattle is the product of ovulation rate, conception rate and embryo survival; twinning and ovulation rate in cattle show a strong, positive genetic correlation (0.75 to 0.9, Van Vleck et al., Reference Van Vleck, Gregory and Echternkamp1991). Heritability estimates for twinning in cattle range from 0.01 in Israeli Holstein heifers unselected for twin births (Weller et al., Reference Weller, Golik, Seroussi, Ron and Ezra2008) to 0.10–0.12 in the USMARC experimental herd selected for twinning (Van Tassell et al., Reference Van Tassell, Van Vleck and Gregory1998). These authors have confirmed that selection for twinning rate is possible despite its low heritability level.
At the molecular level, Vinet et al. (Reference Vinet, Drouilhet, Bodin, Mulsant, Fabre and Phocas2012) presented a detailed review on the genetic basis of multiple births, stating that twinning rate in cattle seems to involve a large number of different genes and different biological pathways. However, strong evidence for the existence of major genes for ovulation rate on BTA 10, where three candidate genes are encoded (SMAD3, SMAD6 and IQCH), has been provided by a recent genome-wide investigation on twinning rate in cattle performed by Kirkpatrick and Morris (Reference Kirkpatrick and Morris2015).
The implicit scope of this study was to investigate the genetic control of twinning in a suckler herd of the Maremmana breed. The first aim was to investigate heritability for twinning, using linear and threshold animal models. A genome-wide analysis was also performed to verify the possible existence of some significant quantitative trait loci for twinning rate. Heritability estimates of twinning rate obtained from pedigree data were then compared with heritability estimates obtained from genotypic data.
Material and methods
Experimental animals and data set
The herd of Maremmana breed used as a case study belongs to the experimental farm of the Animal Production Research Centre in Monterotondo (CREA), near Rome in Italy. The herd has been kept on range all year long in two separate, but genetically connected groups, since the 1930s. In recent years, one group has been maintained purebred to assure replacements, whereas the other group has been crossed with Piedmontese bulls to produce calves with a better beef conformation. The Maremmana is a hardy breed usually kept on range in marginal areas with very poor pastures, and the production system requires seasonal reproductive cycles, calving from February to early May. The breeding bulls – either Maremmana or Piedmontese – are left with the cows roughly from May through July. Calving records, including type (single or twin) were available since 1983. Because of lack of better information, all twins were considered dizygotic. Although dizygotic twins result from multiple ovulations, monozygous twinning results from embryo splitting during early development, for causes unknown so far. In cattle, monozygous twinning is biologically a very rare phenomenon, which occurs in about 2% to 10% of twins of the same sex (Hancock, Reference Hancock1954; Johansson et al., Reference Johansson, Lindhé and Pirchner1974). In any case, since calving occurs outdoor with no human aid, the twinning type is very hard to detect.
The data set included 1260 individuals born from 1963 to 2014, 527 males and 733 females, 402 of them calving at least once from 1983 through 2015. Also included in the data set were 37 sires of dams: three with >50 daughters, eight with 16 to 50 daughters, nine with 6 to 15 daughters and 17 bulls with <6 daughters each.
The number of recorded calvings from 1983 through 2015 was 2130: 2009 single and 121 twin births. Until 2010, herd management did not aim at increasing twinning rate, therefore twins were systematically discarded from reproduction.
During 4 years, from 2012 through 2015, the herd was closely monitored for reproductive events and calf survival.
Statistical analyses
Prediction of breeding values
Prediction of breeding values and estimates of variance components together with their standard errors, were obtained by a single-trait animal model including permanent environmental factor:
where y ijk is the number of calves born at the k th calving of the i th cow; C j fixed effect of the j th year of calving (j=1, 2, … , 33); age ik age in days of the i th cow at its k th calving; F i inbreeding coefficient of the i th cow; b 1, b 2, b 3 regression coefficients; p i permanent environmental effect of the i th cow; a i random additive genetic effect of the i th animal; e ijk random residual effect for the k th calving of the i th cow.
Because twinning is a dichotomous trait with clearly unbalanced frequencies between the two categories, we also fitted (1) by an animal threshold model.
Statistical analyses were performed by the package Multivariate Mixed Models DMU (Derivative-free approach to Multivariate analysis by Restricted Maximum Likelihood), version 6.0, release 5.2 (Madsen and Jensen, Reference Madsen and Jensen2013). The module based on Average Information Restricted Maximum Likelihood (derivative-free approach to multivariate analysis by restricted maximum likelihood average information (DMUAI)) was used to fit linear models, while the module with a Gibbs sampler based on the Reversible Jump Markov Chain (RJMC) algorithm was used to fit threshold models and to estimate variance components.
Because the DMU package does not provide heritability and its standard error for models with more than one random factor, standard error of heritability
$\big( {{{\hat{\sigma }_{a}^{2} } \mathord{\left/ {\vphantom {{\hat{\sigma }_{a}^{2} } {\hat{\sigma }_{p}^{2} }}} \right. \kern-\nulldelimiterspace} {\hat{\sigma }_{p}^{2} }}} \big)$
was estimated following Su et al.’s (2007):
$${\rm SE} \left[ {{{\hat{\sigma }_{a}^{2} } \over {\hat{\sigma }_{p}^{2} }}} \right]{\equals}\sqrt {\left[ {{{\hat{\sigma }_{a}^{2} } \over {\hat{\sigma }_{p}^{2} }}} \right]^{2} {\times}\left[ {{{{\rm Var}(\hat{\sigma }_{a}^{2} )} \over {(\hat{\sigma }_{a}^{2} )^{2} }}{\plus}{{{\rm Var}(\hat{\sigma }_{p}^{2} )} \over {(\hat{\sigma }_{p}^{2} )^{2} }}{\minus}2{{{\rm Cov}(\hat{\sigma }_{a}^{2} , \hat{\sigma }_{p}^{2} )} \over {\hat{\sigma }_{a}^{2} \hat{\sigma }_{p}^{2} }}} \right]} $$
where
$\hat{\sigma }_{a}^{2} $
is the estimated additive genetic variance and
$\hat{\sigma }_{p}^{2} $
estimated phenotypic variance.
Estimates of (co)variance components were obtained from the output of DMUAI or computed from the distributions generated by the RJMC algorithm described above.
Long-term genetic trend was estimated by linear regression of average estimated breeding values (EBVs) per year of birth of animals over time.
Correlation between twinning, fertility and calf mortality
In the Maremmana breed calving may only occur in spring, then calving interval is not the optimal criterion to evaluate fertility. In order to take into account the permanent environmental effect of the cow, a multitrait threshold animal model with the number of newborn per cow in a given year (1 or 2) and the occurrence of calving in the next year (1 or 0) for the same cow was implemented. However, such model did not reach convergence.
Then, in order to verify whether twinning had a negative effect on fertility, a two-trait threshold animal model was used to estimate the genetic correlation between the number of twin calvings during the lifespan of a cow and the number of years each cow stayed empty during lifespan.
The following model was used:
where y ijl for l=1: number of twin calvings of the i th cow during lifespan; y ijl for l=2: number of years the i th cow did not calve during lifespan; C jl fixed effect of number of years from first calving to culling (j=1, 2, … , 17) on the l th trait; F il inbreeding coefficient of the i th cow; b l regression coefficient for the l th trait; a il random additive genetic effect of the i th animal on the l th trait; e ijl random residual effect for the i th cow and the l th trait.
To further investigate the phenotypic correlation between twinning and fertility, the χ 2 test of independence was used to estimate the probability of missing calving in the year following a twin or a single calving.
As regards to calf mortality, from 2012 through 2015 a χ 2 test of independence was used to analyze the relation between twinning and calf mortality, overall and within group, purebred v. Piedmontese crossbred.
Genome-wide association study
Genotyping was performed on the females living on farm in 2012 (119 head) using the BovineSNP54 BeadChip (Illumina, San Diego, CA, USA) which contains 54 524 single nucleotide polymorphisms (SNPs). Although acknowledging the limited sample size, a genome-wide association analysis of twinning ability was performed by the GenABEL package in R (Aulchenko et al., Reference Aulchenko, Ripke, Isaacs and van Duijn2007) using a three step GRAMMAR-GC approach (genome-wide association using mixed model and regression – genomic control). In fact, this procedure estimates the genomic kinship matrix from genomic marker data, so allowing to compare heritability values obtained through pedigree with genomic heritability. On the other hand, considering the scant literature on molecular markers for twinning in cattle, the detection of genomic regions possibly associated with this trait might add new insight into the mechanisms underlying the trait. In order to perform genomic analyses by the GenABEL software, it was necessary to define a unique phenotype per each genotyped cow. EBVs could be used, but their use in association studies with DNA markers could result in higher level of false positives (Ekine et al., Reference Ekine, Rowe, Bishop and de Koning2014). Therefore uncorrelated random effects of each cow obtained by the following model with no consideration of relationships among animals (assuming V=
I
$\sigma _{a}^{2} $
) were used in the genome-wide association study (GWAS).
where y ijk is the number of calves born at the k th calving of i th cow; C j fixed effect of the j th year of calving (j=1, 2, … , 33); age ik age in days of the i th cow at its k th calving; F i inbreeding coefficient of the i th cow; b 1, b 2, b 3 regression coefficients; a i uncorrelated random effect of the i th cow; e ijk random residual effect for the k th calving of the i th cow.
This procedure had been previously applied for the estimation of the effects of molecular markers in presence of repeated measures for each individual (Crisà et al., Reference Crisà, Marchitelli, Pariset, Contarini, Signorelli, Napolitano, Catillo, Valentini and Moioli2010; Moioli et al., Reference Moioli, Scatà, Steri, Napolitano and Catillo2013).
Markers were excluded from the analysis if minor allele frequency (MAF) was <1% and call rate <95%. Call rate for samples was also computed and cows with more than 10% missing data were removed.
Classical Multidimensional scaling (MDS) was performed within the R statistical environment using the GenABEL package (Aulchenko et al., Reference Aulchenko, Ripke, Isaacs and van Duijn2007) to explore population structure and verify the genetic homogeneity of the data set before carrying out the GWAS.
Results
Twinning and calf mortality
From 2012 through 2015, twinning rate grew from 5.63% to 12.77% (Table 1).
Twinning rate from 2012 to 2015

During the 4 years, total mortality was 8.93%, with significant differences between single and twins (5.23% v. 26.57%; P<0.0001) as well as between crosses and purebred Maremmana (4.40% v. 12.77%; P=0.006) (Table 2).
Calf mortality in single and twin births, in the Maremmana (Mar) purebred v. Piedmontese crossbred (Pied)
Total mortality in purebred single calving (8.44%) was significantly higher than in crossbred single calving (1.50%; P=0.009), whereas mortality in purebred twin calving (32.35%) was not significantly higher than crossbred twin calving (19.23).
Stillbirth in single purebred calving (3.90%) was not significantly higher than in single crossbred calving (0.75%), whereas stillbirths were very similar in twin purebred calving (14.71%) and in twin crossbred calving (15.38%).
Mortality from birth to weaning in single purebred calving (4.55%) was significantly higher than in single crossbred calving (0.75%; P=0.052). Mortality from birth to weaning in twin purebred calving (17.65%) was not significantly higher than twin crossbred calving (3.85%). However, considering together single and twin calving, mortality from birth to weaning was significantly higher in purebred (7.60%) than in crosses (1.26%; P=0.01).
Only two difficult calvings were recorded, both of them in single calving due to breech presentation. One twin calf was abandoned by its dam, it was taken in the barn and artificially fed until weaning.
Breeding values and genetic parameters for twinning
Table 3 reports the genetic parameters estimated by the linear and the threshold models; heritability by the threshold model (0.062±0.093) was higher than by the linear model (0.014±0.018). Descriptive statistics on EBVs for twinning calculated with both models are reported in Table 4.
Genetic parameters for twinning (SE in parentheses) estimated by Repeatability Linear Model (R_LM) and Repeatability Threshold Model (R_ThM)

$\sigma _{a}^{2} $
=additive genetic variance;
$\sigma _{{pe}}^{2} $
=permanent environmental variance;
$\sigma _{e}^{2} $
=error variances; r=repeatability:
$${{\big( {\sigma _{{pe}}^{2} {\plus}\sigma _{a}^{2} } \big)} \mathord{\left/ {\vphantom {{\left( {\sigma _{{pe}}^{2} {\plus}\sigma _{a}^{2} } \right)} {\left( {\sigma _{{pe}}^{2} {\plus}\sigma _{a}^{2} {\plus}\sigma _{e}^{2} } \right)}}} \right. \kern-\nulldelimiterspace} {\big( {\sigma _{{pe}}^{2} {\plus}\sigma _{a}^{2} {\plus}\sigma _{e}^{2} } \big)}}$$
; h
2=heritability:
$${{\left( {\sigma _{a}^{2} } \right)} \mathord{\left/ {\vphantom {{\left( {\sigma _{a}^{2} } \right)} {\left( {\sigma _{{pe}}^{2} {\plus}\sigma _{a}^{2} {\plus}\sigma _{e}^{2} } \right)}}} \right. \kern-\nulldelimiterspace} {\big( {\sigma _{{pe}}^{2} {\plus}\sigma _{a}^{2} {\plus}\sigma _{e}^{2} } \big)}}$$
.
EBVs estimated by Repeatability Linear Model (R_LM) and Repeatability Threshold Model (R_ThM)

Figure 1 reports average EBVs per year of birth predicted by the threshold model. Genetic trend for twinning over a 50-year span (1963 to 2014) was –0.0001, null for all practical purposes. This is compatible with the lack of interest towards increasing twinning rate, together with the low frequency of the event.
Yearly averages of EBVs for twinning rate from 1963 to 2014.
Correlation between twinning and fertility
Out of 1945 single births available in the data set, 1590 were followed by a calving in the next year (81.7%), while 96 out of 120 twin births were followed by a calving in the next year (80.0%); χ 2 was evidently not significant.
The correlation between additive genetic and permanent environmental effects for twinning rate and number of years each cow did not produce calves was −0.27±0.03.
Genome-wide association study
After data editing, 409 markers having SNP call rate <0.95 and 12 584 markers having MAF <0.01 were excluded, while 6 cows had to be removed after the implementation of the individual call rate of 0.90. No outliers were identified by MDS (Figure 2) consequently the data set passing quality controls and used in the association analysis contained 41 638 SNPs and 113 samples. Mean heterozygosity of the sample was 0.363±0.037 for individuals, and 0.337±0.148 for SNPs. Genomic heritability, estimated with the ‘ibs’ function of GenABEL (option weight=‘freq’), was 0.29±0.021. Evidence of population stratification was estimated by the genomic inflation factor λ=1.02 for a basic χ 2 statistics test, and was completely corrected by the GRAMMAR-GC methodology that yielded λ=1. The Q-Q plot (Figure 3) shows no systematic deviation from the diagonal indicating that the data set was corrected for population stratification. Associated SNPs deviate from the diagonal at the upper-right end of the plot.
Multidimensional scale plot of the entire sample set of animals included in the association analysis.
Q-Q plot of (GRAMMAR-GC) P-values/against expected P-values.
Manhattan plots displaying the GWAS results with respect to their genomic position are shown in Figure 4. The threshold for significance was set according to The Wellcome Trust Case Control Consortium (2007) and to Panagiotu and Ioannidis (Reference Panagiotu and Ioannidis2012). One significant SNP (Hapmap22923-BTA-129564) was identified on chromosome 24 at position 40020911 bp with P-value of 6.47 e −06 and MAF of 1.32%; the SNP explained 0.034% of the phenotypic variance of twinning. Because of the limited sample size, the search for putative candidate genes in linkage disequilibrium (LD) with the significant SNP was restricted to a genomic region spanning ±300 kb around the marker, even though cattle breeds are characterized by large LD blocks, likely caused by evolutionary and selective forces (Odani et al., Reference Odani, Narita, Watanabe, Yokouchi, Sugimoto, Fujita, Oguni, Matsumoto and Sasaki2006; Khatkar et al., Reference Khatkar, Zenger, Hobbas, Hawken, Cavanagh, Barris, McClintock, McClintock, Thomson, Tier, Nicholas and Raadsma2007). Two annotated genes were detected on chromosome 24, in proximity of Hapmap22923-BTA-129564: ARHGAP28 gene, 208 kb downstream the marker and TMEM200C gene 282 kb upstream.
Manhattan plot displaying the results (−log10 of P-value) of the genome-wide scan using the GRAMMAR-GC method with respect to the genomic position. The horizontal lines separate single nucleotide polymorphisms (SNPs) with P-values <0.001 and <0.01, respectively.
Discussion
Today, selection objectives in animal farming pay a great deal of attention to robustness, longevity, welfare and, more generally, sustainability of production. Negative factors associated with twinning, as reduced calf survival, increased incidence of dystocia, and longer intervals between conceptions (Kirkpatrick, Reference Kirkpatrick2002) have generally prevented any increase in twinning ratio, mainly in dairy production systems. The situation may be different for beef cows where the higher production of twin calves may represent an opportunity to increase profitability within a sustainable livestock system (Echternkamp and Gregory, Reference Echternkamp and Gregory1999). Van Vleck and Gregory (Reference Van Vleck and Gregory1996) reported an increase in twinning rate from 3.4% in 1982 to 28.5% in 1993 after selecting for twinning. This study on the Maremmana breed showed that this selection goal is worth to be pursued, even if a moderate negative genetic correlation was detected between twinning rate and the number of years cows did not produce calves; despite of the referred association between twinning and fertility, to our knowledge no estimate of genetic correlation between the two traits has been reported in literature. At the phenotypic level, the hypothesis that a cow will be more likely empty in the year following a twin calving, was here rejected. This phenotypic result is highly informative because it comes from data collected on cows with very high longevity (replacement rate <10%).
Higher mortality of twin calves was the only detected adverse effect; however, the number of extra dead calves was lower than the increase in newborns and the incidence of calf mortality was lower in crossbred. At any rate, improvements in herd management may further reduce calf mortality. In the present study, twin births increased from 2012 through 2015 (Table 1) likely for environmental reasons since selection for twinning had not started yet.
Twinning is a dichotomous trait with clearly unbalanced frequencies between the two categories, and in such cases threshold models appear to be the method of choice for data analysis (Moorad and Promislow, Reference Moorad and Promislow2011). Heritability of dichotomous traits estimated by threshold models are higher than from linear models (Ron et al., Reference Ron, Ezra and Weller1990; Ghavi Hossein-Zadeh et al., Reference Ghavi Hossein-Zadeh, Nejati-Javaremi, Miraei-Ashtiani and Kohram2009) and this was fully confirmed in this study. Heritability for twinning in this study was generally low, confirming heritability values estimated worldwide. Van Vleck and Gregory (Reference Van Vleck and Gregory1996) and Gregory et al. (Reference Gregory, Bennett, van Vleck, Echternkamp and Cundiff1997) reported heritability values of 0.02 to 0.08, Syrstad (Reference Syrstad1984) found in Norwegian dairy cattle values ranging from 0.006 in the first parity up to ~0.04 in parities 3 to 5. Finally, Ron et al. (Reference Ron, Ezra and Weller1990) reported heritability values of 0.022 from a sire linear model and 0.101 from a threshold sire model on 124 553 calving records in Israeli Holsteins. Later studies (Ghavi Hossein-Zadeh et al., Reference Ghavi Hossein-Zadeh, Nejati-Javaremi, Miraei-Ashtiani and Kohram2009) on larger number of Holstein (1 163 594 calving events) reported heritability values for parities 1, 2 and 3 of 0.018, 0.014 and 0.011, respectively, by an animal linear model and 0.074, 0.067 and 0.054 by an animal threshold model. The estimates by those authors were very similar to the results obtained in this study from a much smaller sample. Weller et al. (Reference Weller, Golik, Seroussi, Ron and Ezra2008) found heritability values ranging from 0.01 to 0.03 using an animal linear model on Israeli Holstein heifers unselected for twin births.
Genomic heritability calculated in this study by the genotype kinship matrix for cow twinning ability was much higher (0.29) than any heritability estimate through pedigree registrations. According to simulation studies (Zaitlen et al., Reference Zaitlen, Pasaniuc, Sankararaman, Bhatia, Zhang, Gusev, Young, Tandon, Pollack, Vilhjálmsson, Assimes, Berndt, Blot, Chanock, Franceschini, Goodman, He, Hennis, Hsing, Ingles, Isaacs, Kittles, Klein, Lange, Nemesure, Patterson, Reich, Rybicki, Stanford, Stevens, Strom, Whitsel, Witte, Xu, Haiman, Wilson, Kooperberg, Stram, Reiner, Tang and Price2014) heritability estimated under the assumption of the infinitesimal model tends to be higher than the genomic heritability. However, in this experiment, heritability estimated for the entire data set was much smaller than heritability estimated for the subset of genotyped cows. This result may depend on: (i) the dichotomous nature of the trait; (ii) the strong prevalence of single births; (iii) the genetic control of the trait, possibly based on a limited number of genes. On the other hand, the genomic kinship should better represent the true relations between two individuals, because it is based on the alleles that each individual actually carries and not on the probability of receiving those alleles by descent. Although genomic heritability was calculated on 113 individuals only (i.e. the genotyped cows), it provides further clues on the possibility of selecting for this trait. Because the study dealt with the feasibility of on farm selection of a local breed, sample numbers were limited, then the association between SNPs and twinning rate should not be considered the predominant result of this study. This also because the GWAS was performed on a ‘corrected phenotype’ including all effects permanently tied to the cow. Therefore only the DNA region around the most significant marker (Hapmap22923-BTA-129564) was explored, and two potential candidate genes for twinning were detected: ARHGAP28, coding for the Rho GTPase activating protein 28, which regulates actin stress fibre assembly; TMEM200C, coding for the transmembrane protein 200C, which is mostly involved in the communication of the cell with the external environment. However, only for ARHGAP28 a direct involvement in the ovulation process might be hypothesized. This gene was identified by Catalano et al. (Reference Catalano, Critchley, Heikinheimo, Baird, Hapangama, Sherwin, Charnock-Jones, Smith and Sharkey2007) among the progesterone-regulated genes in humans and mice, because it showed a −3.1-fold change in gene expression from women treated with antiprogestin as compared with control.
Conclusion
As reported by Komisarek and Dorynek (Reference Komisarek and Dorynek2002), the production of twin calves offers the opportunity to increase cow economic efficiency while entailing some negative repercussions. These authors estimated that, even taking into consideration higher labour and veterinary costs, the increase in efficiency of beef production due to twinning can reach 24%. The present study confirmed the possibility to control twinning rate in cattle, because this parameter showed some genomic heritability; moreover, the potential of twinning to increase the productivity of local breeds like the Maremmana, kept in production systems where other reproductive traits could hardly be improved, was here described. Unfavourable effects associated to twin births should be taken into account when selecting for twinning rate: in order to reduce death rate of calves before weaning, consideration should be given to mothering capability and management improvement. However, deaths more frequently occurred in purebred calves, and this result supports the existence of an effect of the calf genotype. During this study, dystocia in the Maremmana breed was not a problem, even in crosses with the Piedmontese beef breed. Calf abandoning was rarely registered during the first 10 days after calving, and it could be easily controlled by careful management. Although heritability computed from pedigree data was very low, genomic heritability was higher and this might be a promising hint (signal) that twinning rate could be increased even in small cattle populations, also explaining the relevant selection response obtained in different experiments around the world. Association studies on larger number of animals could allow for deeper investigations in the genomic control of twinning. Finally, in systems where it is not a desirable event, selection against twinning should also be possible.
Acknowledgements
This research was supported by the governmental research programme 19735/7303/2012 ‘REDDBOV’ provided by the Italian Ministry of Agriculture (Rome, Italy).


