Ear size variation traits analyses
The ear size among different cattle individuals and breeds were obviously different. As a general observation, the ear size of Bos indicus in Indo-Pak region was rather bigger than the Bos taurus in the world, while hybrid cattle breeds depicted diversified ear size (Fig. 1a). According to this observational phenomenon, the ear size data was collected and calculated from 475 adult female cattle include 6 representative breeds(i.e. Burmese, Brahman, Simmental, Yunling, Wenshan and Dabieshan cattle) using pixel method (29). The overall ear size data from 87.81 to 330.85 cm2. Bos indicus (Brahman and Burmese cattle) have the largest mean values, which were 226.75 cm2 and 254.06 cm2, respectively. In contrast, Bos taurus (Simmental) has the small mean value, which was 134.62 cm2. The size of hybrid breed (Yunling) and Chinese indicine (Wenshan and Dabieshan cattle) were relatively scattered, with an average values 164.21, 186.31, and 184.37, respectively (Table 1). The ear size distribution of all individuals showed a unimodal distribution according the histogram and density plot (Fig. 1b). And it can be well distinguished in three distinct intervals: big ear (> 230 cm2), middle ear (< 230 cm2 & >148 cm2) and small ear (< 148 cm2) with 25% and 75% quantiles of the all individuals’ boxplot analysis (Fig. 1c). Boxplot and F statistic was used to analyze the variance in the ear size among six cattle breeds. The results showed that Bos indicus depicted largest ear size, while Bos taurus represented the smallest ear size, while the hybrid (Yunling) and the Chinese indicine (Wenshan and Dabieshan cattle) ranged in the middle (Fig. 1c). Moreover, the F statistic results also suggested that the ear size among cattle breeds were significantly different, supporting the above statement.
Table 1
Descriptive statistics of ear size in cattle breeds
Breeds | Origin | Number | Maxmum | Minmum | Mean | CV |
Simmental | Bos Taurus | 65 | 166.26 | 92.67 | 134.62 | 10.53% |
Brahman | Bos indicus | 39 | 330.85 | 140.86 | 226.75 | 15.68% |
Burmese | Bos indicus | 52 | 313.53 | 212.46 | 254.06 | 8.50% |
Yunling | Bos indicus × Bos Taurus | 119 | 268.71 | 87.81 | 164.21 | 18.88% |
Wenshan | Chinese indicine | 100 | 315.87 | 95.47 | 186.31 | 22.28% |
Dabie mountain | Chinese indicine | 100 | 249.86 | 120.63 | 184.37 | 13.75% |
aCoefficients of variation |
Two genome-wide association studies for ear size traits
Based on 13,057,965 autosomal SNPs derived from the 158 published resequencing cattle data (18), a genome-wide efficient mixed-model association (GEMMA) (23) was used in the primary genome-wide association studies to identify the significant loci. Figure 2a showed the Manhattan plot for the GWAS. The red and green horizontal lines represent the Bonferroni-adjusted genome-wide significant and suggestive threshold. A total of 3 significant SNPs and 40 suggestive SNPs showed genome-wide associations with ear size (Table S3). However, most of the significant potential SNPs were located on BTA 6 (36.79 ~ 38.90 Mb) (Fig. 2b). The significant locus (P = 1.29 × 10− 8) was located in LOC112447053, and near MEPE and IBSP. (Fig. 2c). The Quantile-Quantile plot (QQ-plot) in Fig. 2d showed the observed and expected P-values of the GWAS for ear size. The dashed line represents the distribution of the SNPs under the null hypothesis that there is no association of SNPs with the trait of interest. The strong deviation of the observed from the expected P-values for QQ-plots indicate that there were more SNPs significantly associated with all of the ear size trait than would be expected by chance.
Since, body size traits and ear size data could have a shared underlying genetic basis, the correlation was measured and analyzed for each trait in Figure S1. Pearson correlation coefficient was used to assess the linear relationships of ear size and body size traits. The result showed that the body height, cross high, head length and head width had a weak-to-moderate positive correlation with the ear size trait. The multiple linear regression model was analyzed using these related traits as covariates for secondary GWAS analysis using PLINK (24). The autosomal SNP scan for ear size revealed associated markers as showed in Figure S2. A total of 293 SNPs were observed on the potential region and most of them were found on BTA6. The most significant SNP (P = 5.74 × 10− 11) was located on the intron of IBSP. The two GWAS results strongly suggested that IBSP had a strong correlation with ear size and might be a key candidate gene influencing cattle ear size.
Genetic differentiation of the mutation in IBSP among different ear size groups
To clarify whether ear size was selected in different cattle breeds, pairwise fixation index (Fst)(30) was used to measure the genetic differentiation between big (Brahman, Burmese) and small ear group (Simmental, Yanbian cattle). By annotating the significant regions (Fst > 0.66, empirical Pvalue < 0.005), IBSP was found in the most significant region on BTA6 (position: 36,840,001 ~ 36,890,000, Fst = 0.75) (Fig. 3a). To further investigate the differentiation of these mutation loci across diverse cattle breeds, 14 cattle breeds were used for nucleotide diversity (θπ) and SweepD analysis (Figure S3, S4) (31, 32). The results showed that the SweepD value on these loci were observed to be near zero in Bos taurus (Angus, Kazakh, Hanwoo, Hereford, Holstein, Mishima, Simental, and Yanbian cattle) and over 200 in South Asian indicine (Brahman, Burmese), which means that IBSP may be selected in indicine breeds (Fig. 3b). According to the functional annotation, a missense mutation was found on the seventh exon of IBSP (Threonine-250→Isoleucine, T250I) (Fig. 3c). To further evaluate the functional impact of the variants, we aligned the mutant IBSP protein with its ortholog proteins in bovidae (Fig. 3d) and other diverse vertebrates (Figure S5). The comparison reveals that T250I was a quite conserved amino acid mutation, which is invariant among all the other animals we examined except Bos Taurus. We calculated the linkage disequilibrium (LD) values of the SNPs, which were shared in the region of IBSP. We observed strong linkage in this region (Figure S6). At the same time, prediction based 3D structure of protein showed that the mutation site of IBSP can change its protein structure (Fig. 3e). All those results imply that T250I is most likely causal mutation for the IBSP sweep in cattle ear size.
Allele frequency of IBSP mutation among different cattle breeds indicated the origin
In our study, we genotyped the DNA sequences of IBSP with 394 cattle representing Bos indicus, Bos taurus and hybrid cattle to investigate the genotype frequency of the mutation loci across diverse cattle breeds. Two genotypes (A and G) of the IBSP mutation loci were found (Fig. 4a, 4b). The A allele of IBSP mutation occurred at high frequency in Burmese (0.06), Brahman (0.30) and Lingnan (0.40), respectively. In Wenshan, Shigatse, Weining, Hainan, Dianzhong, Jinjiang, Ji'an, Wannan and Luxi cattle, the A allele of the two mutation showed a moderate frequency close to 0.5. In contrast, the A allele was almost equal to one in Yanbian, Kazakh, Simmental and other Bos taurus (Fig. 4a, Table S4). At the same time, we searched the mutation frequency of IBSP among different cattle breeds in the world from the Bovine Genome Variation Database and Selective Signatures (BGVD, http://animal.nwsuaf.edu.cn/code/index.php/BosVar) (33). The database contained 24 South Asian indicine cattle, 19 Chinese indicine, 37 East Asian taurine, 38 European taurine, 19 Eurasian taurine and 10 Africa taurine cattle. The results showed that the A allele frequency of IBSP was highest in South Asian indicine (0.10), and the frequency was 0.47 in Chinese indicine, respectively. However, the A frequency of IBSP was near or equal to one in East Asian taurine, European taurine and Eurasian taurine (Fig. 4b). The allele distribution indicates that the A allele of IBSP may have originated in Bos taurus. Due to the infiltration of South Asian indicine, the Chinese indicine caused the reduce in the frequency of A allele. We also found the G allele in African taurine at low frequency, which may be due to the large introduction of South Asian indicine after the African rinderpest in 1890 (34).