Genomic diversity within and among the breeds
Twenty of the 132 animals were excluded due to low genotyping (MIND > 0.1), and the average genotyping rate for the remaining 112 animals was 0.99. The final data on 112 cattle samples belonging to Sahiwal (13), Tharparkar (17), Gir (15), Ongole (17), Hariana (18), Kangayam (16) and Vechur (16) breeds were achieved after quality control measures (Table 1). A sample size of ≥ 13 per breed was adequate for the diversity analyses, which was in consonance with other studies18,19. Out of 735293 autosomal SNPs genotyped in these cattle, 60% (438,176) were in Hardy Weinberg equilibrium and revealed higher degree of polymorphism (MAF > 0.05). It was also stated that a sample size of 4-620 and polymorphic SNP filtration19,21,17 could mitigate the effect of ascertainment bias if the number of markers is sufficiently large enough, as in the current study. The SNPs in strong LD (R = 0.5) were also filtered out to minimize the bias in estimating the genomic diversity. We were left with 165,021 informative SNPs after pruning for further genome diversity analyses. Minor allele frequency varied between breeds, ranging from 0.23 (Kangayam) to 0.26 (Vechur), with an average observed heterozygosity of 0.35 in all samples examined (Table 1). The similar estimates of MAF were also reported in Nellore, Holstein, Iranian, Ethiopian and South African cattle22–26, which varied from 0.21 to 0.25 in Bos indicus and Bos Taurus cattle breeds. However, higher/lower estimates were also reported in other cattle breeds27–29. The reported estimates of observed heterozygosity in several other zebu cattle25,30,31 were in agreement with the present estimates. However, higher heterozygosity was reported in Hanwoo, Rwanda and other Taurine cattle breeds25,30,29. The lower genetic variability estimates in Bos indicus relative to Bos taurus were consistent with earlier studies32,33,29. The lower genetic diversity in Bos Indicus may be due to ascertainment bias of mainly taurine derived SNPs in SNP chip. Out of Bovine HD and 54K SNP chips of Illumina, only 40–50% SNPs were found to be informative for genetic diversity of Zebu cattle breeds of India16 but, it was 90% in Bos taurus after quality control14.
Table 1
Number of animals, mean of expected and observed heterozygosity (He, Ho), minor allele frequency (MAF) and coefficient of inbreeding (FIS).
Breed
|
n
|
He
|
Ho
|
MAF
|
FIS (p > 0.40)
|
Sahiwal
|
13
|
0.34
|
0.35
|
0.25
|
-0.009
|
Tharparkar
|
17
|
0.33
|
0.34
|
0.24
|
-0.003
|
Gir
|
15
|
0.32
|
0.33
|
0.24
|
0.010
|
Ongole
|
17
|
0.33
|
0.34
|
0.24
|
0.006
|
Hariana
|
18
|
0.33
|
0.34
|
0.25
|
-0.001
|
Kangayam
|
16
|
0.30
|
0.33
|
0.23
|
-0.062
|
Vechur
|
16
|
0.34
|
0.35
|
0.26
|
-0.009
|
overall
|
112
|
0.36
|
0.35
|
0.27
|
-0.009
|
The genome diversity of seven cattle breeds viz. Tharparkar, Sahiwal, Gir, Vechur, Ongole, Kangayam and Hariana, was determined first using all filtered SNPs (4 38,176) and then using the 1993 most informative SNPs having FST value ≥ 0.20. The genetic differentiation (FST) of the breeds based on all SNPs was just 0.05 and based on 1993 most differentiating loci with an FST value of ≥ 0.20 was 0.23. Hence, out of 777K, 1993 loci (0.25% of the total) were identified as having greater power to distinguish these cattle populations, and these loci were further assessed for diversity and structure analysis. The overall estimates of observed heterozygosity (Ho), gene diversity within population (Hs), gene diversity in total population (Ht), gene diversity among breeds (Dst), and Dest, a measure of population differentiation34, across these loci were 0.28 0.27, 0.35, 0.08 and 0.13, respectively.
Manhattan graph (Fig. 1) depicts the distribution of FST values across the chromosomes. The spread of FST showed that there are few loci (16 markers) having the higher degree of genetic differentiation (FST ≥ 0.40). The inbreeding coefficient (identical by state) was zero, suggesting that the cattle from which the samples were taken were randomly mated. The genetic differentiation power of those informative SNP loci (n = 1,993) in cattle breeds under study ranged from 0.20 to 0.51, with an average value of 0.23, indicating that these loci account for 23% of the genetic variation between the breeds. Table 2 and Fig. 2 showed the pair-wise estimates of F-statistics (FST) and Nei's genetic distance35. Both estimates revealed higher genetic differentiation between Kangayam and rest of the cattle breeds (FST: 0.08 to 0.10), followed by Ongole and rest of the breeds (FST: 0.05 to 0.07). Dairy cattle breeds viz., Tharparkar, Sahiwal and Gir showed the least genetic differentiation (FST: 0.04 to 0.06). A moderate genetic differentiation36 where FST ranged from 0.05 to 0.15 was also recorded in Bos taurus dairy and beef cattle with an FST value of 0.08 across all SNPs14, as well as in African cattle breeds (FST : 0.04 to 0.08). Many SNPs with higher genetic differentiation power (FST > 0.5) were also reported in Bos taurus dairy and beef cattle14. The degree of genetic differentiation among Indian Zebu cattle was comparable to that of African Zebu cattle30, but lower than that of Bos taurus and Bos indicus30,33, which was anticipated given their historical divergence.
Table 2
The pair-wise Nei's FST (above diagonal) and genetic distance (below diagonal) among the cattle breeds
Breed
|
TP
|
SW
|
GR
|
VC
|
OG
|
KG
|
HR
|
TP
|
-
|
0.044
|
0.059
|
0.066
|
0.064
|
0.094
|
0.049
|
SW
|
0.028
|
-
|
0.051
|
0.051
|
0.052
|
0.081
|
0.037
|
GR
|
0.032
|
0.032
|
-
|
0.072
|
0.072
|
0.101
|
0.056
|
VC
|
0.034
|
0.031
|
0.037
|
-
|
0.069
|
0.096
|
0.061
|
OG
|
0.033
|
0.031
|
0.037
|
0.035
|
-
|
0.092
|
0.058
|
KG
|
0.044
|
0.042
|
0.048
|
0.046
|
0.043
|
-
|
0.088
|
HR
|
0.027
|
0.026
|
0.031
|
0.032
|
0.031
|
0.043
|
-
|
SW-Sahiwal, TP-Tharparkar, GR-Gir, VC-Vechur, OG-Ongole, KG-Kangayam, HR-Hariana |
Effective Population Size
Ancestral and recent effective population sizes (Ne) for seven Indian cattle breeds are presented in Fig. 3. Estimated Ne showed a downward trend in recent generations across the populations. The most rapidly declining recent Ne was found in the KG and GR, whereas HR and TP showed a slowly declining trend. The estimated Ne at 13 generations ago in seven Indian native cattle breeds ranged from 48 to 74 (Supplementary Table S1).
Over the last 999 years, the Ne has displayed a decreasing pattern across Indian breeds, with a steeper slope since about 200 generations ago. Lower Ne had been found in KG and SW with estimated values of 48 and 51, respectively, 13 generations ago (Ne13), due to the intensive selection pressure or artificial insemination used for developing these breeds. In Kangayam cattle, we have observed maximum autozygosity attributable to both recent and ancient inbreeding17.
Recently, a Bovine HD-SNPs array based screening of Chinese native cattle populations revealed a similar trend, with Ne13 values ranging from 85 to 13237. The genotyping of Italian local beef breeds (Calvana-CAL, Mucca Pisana-MUP, Pontremolese-PON) and Italian Limousin (LIM) using the GeneSeek GGP-LDv4 33k SNP chip containing 30,111 SNPs showed an average estimated historical effective population size (Ne13) of 45–310 (CAL-79, MUP-65, PON-45 and LIM-310)38.The application of LD-based Ne estimation in developing countries for local breeds without pedigree information could offer new perspectives for the assessment of the actual gene pool available and the respective decision-making in conservation and management.
Genetic Structuring And Classification Of The Cattle Breeds
The principal component analysis (PCA) based on genomic relationship (IBS-Identity by state) matrix using 165,021 LD pruned autosomal SNPs was undertaken to assess breed composition of the animals. The first and second principal components accounted for 29% and 12 % of the total variation, respectively. Figure 4A, B presented the first, second and third principal components where Vechur, Kangayam and Ongole clustered separately from other breeds under investigation, and rest of the breeds (Gir, Sahiwal, Tharparkar and Hariana) grouped together. The structure and PCA both revealed that dairy breeds clustered apart from the dual (Ongole) except Hariana, and draft breeds. PCA based on SNP data clearly separated the breeds by their utility and size, which was in consonance with analysis of molecular variation (Table 3). Kangayam, a draft breed and Vechur, a small statured breed were quite distinct from the dairy and dual breeds and was in agreement with structuring of these cattle based on runs of homozygosity (ROH)17. Moreover, the structure and genomic relationship among these cattle was also studied using a small number of highly differentiated loci (180, 50, 10 markers) and it was interesting to note that even 10 markers also revealed the same level of differentiation as with large sized battery of markers. Therefore, the structure and relationship between the cattle breeds were not significantly influenced by ascertainment bias, which was consistent with a previous study in sheep39 and cattle17.
Table 3
Hierarchical analysis of variance based on pairwise differences*
Source of variance Variance
|
d.f.
|
Sum of squares
|
variation components
|
Percentage of variation
|
p value
|
Among groups
|
3
|
314971.297
|
486.39137
|
1.64
|
0.029
|
Among populations within groups
|
3
|
234688.652
|
1613.69233
|
5.44
|
0.00
|
Among individuals
within populations
|
105
|
2870271.011
|
-243.80843
|
-0.82
|
0.645
|
Within individuals
|
112
|
3116235.500
|
27823.53125
|
93.75
|
0.001
|
*Significance tests based on 1023 permutations |
The Hierarchical F-statistics computed across different groups indicated significant differences in milk production categories (High: dairy breeds; Moderate: Dual breeds; Low: Draft and small statured breeds) and hence differentiated dairy (Sahiwal, Gir and Tharparkar), dual (Ongole, Hariana), and draft (Kangayam) and small sized cattle (Vechur) (Table 4). PCA based on genotype displayed similar categorization of these cattle17. The analysis of molecular variance also revealed significant (p < 0.05) differences in the proportion of variation (1.64%) due to their functional characteristics viz. dairy, dual, draft and small size. However, there were significant differences (5.44%) between the breeds having a particular characteristic like dairy and dual (Table 3). ROH regions observed in these cattle breeds were also able to differentiate dairy and draft breeds as well as small stature cattle17.
Table 4
Hierarchical F-statistics computed over geographical distribution, body size and level of milk production of the 7 cattle breeds
Hierarchical Level
|
Geographical Region
|
Body Size
|
Level of dairy performance
|
Individual
|
Geographical Region
|
0.00
|
0.19
|
0.44
|
0.43
|
Body Size
|
0.00
|
0.00
|
0.30**
|
0.29
|
Level of dairy performance (high, medium and low)
|
0.00
|
0.00
|
0.00
|
-0.02
|
**p = 0.001 |
Admixture Analyses Of The Cattle Breeds
The population structure of these cattle was studied using admixture model based clustering implemented in Structure40 to partition genome of each animal into a predefined number of clusters. The breeds were clearly grouped into draft and dairy/dual breeds (Fig. 5) for pre-defined K = 2, indicating shared ancestry of dairy and dual purpose breeds (Gir, Tharparkar, Sahiwal, Vechur, Ongole, Hariana) very similar to sharing of paternal lineages41. More than 90% of genome of dairy /dual cattle except Ongole (86%) clustered together and 90% of the genome of draft cattle Kangayam grouped into separate cluster. At K = 3, Sahiwal, Tharparkar, Gir, Vechur and Hariana (> 72%) clustered together. Kangayam and Ongole clustered into their own group (89%). All major dairy breeds (Tharparkar, Gir and Sahiwal) share > 80% of their genome. However, at K = 4, Vechur, was clustered apart from other breeds and at K = 5, Gir also clustered separately. At K = 7, most of animals (> 76%) of the different populations under study except Vechur clustered into their own group of animals called breed (Fig. 5). Vechur population displayed highest rate of admixture (Table 5), consistent with its crossing with other breeds and the same was also revealed based on its ROH analysis17. Ongole, Kangayam and Hariana shared comparatively less of their genome (≤ 15%) with other breeds. The Fig. 6 revealed the increasing mean value of log likelihood across the inferred clusters and indicated the all seven breeds evolved from their independent ancestry. However, the structure at different values of K indicated intermixing of these breeds during the recent past. The admixture analysis (k = 2) revealed the structure of these cattle except Vechur similar to that unraveled by principal component analysis.
Table 5
Proportion of membership of each pre-defined population in each of the 7 clusters
Population
|
Given Inferred Clusters
|
Population size
|
|
1
OG
|
2
HR
|
3
KG
|
4
GR
|
5
TH
|
6
VC
|
7
SW
|
|
TH
|
0.013
|
0.107
|
0.014
|
0.034
|
0.766
|
0.016
|
0.049
|
17
|
SW
|
0.014
|
0.095
|
0.014
|
0.035
|
0.047
|
0.026
|
0.768
|
13
|
GR
|
0.014
|
0.073
|
0.014
|
0.764
|
0.064
|
0.011
|
0.06
|
15
|
VC
|
0.033
|
0.07
|
0.052
|
0.026
|
0.028
|
0.733
|
0.058
|
16
|
OG
|
0.854
|
0.082
|
0.01
|
0.002
|
0.015
|
0.011
|
0.027
|
17
|
KG
|
0.026
|
0.024
|
0.889
|
0.002
|
0.027
|
0.01
|
0.022
|
16
|
HR
|
0.014
|
0.89
|
0.005
|
0.02
|
0.026
|
0.01
|
0.035
|
18
|
SW-Sahiwal, TP-Tharparkar, GR-Gir, VC-Vechur, OG-Ongole, KG-Kangayam, HR-Hariana |
Selection Signatures Among The Breeds
The genomic regions showing most differentiation among breed pairs based on FST value were identified. The most differentiated breed pairs based on their performance were chosen for the identification of selection signatures. Therefore, Kangayam, a draft breed, was used as control breed for the analysis of dairy breeds while, Vechur, a small sized breed for the rest of large sized breed to search for signatures that may be associated with stature. Genome wide significance level was set to 0.001 to represent a selection signatures. Out of the significantly differentiated genomic regions, the top most differentiating signatures with FST ≥0.50 among different breed pairs is presented in Table 6.
Table 6
The selection signatures (FST ≥0.50) value among the different cattle breeds and their annotation with Bovine QTL database
Breed** pair
|
Top FST value
|
Marker
|
Chromosome
|
QTLs*
|
SW/TP
|
0.67
|
BovineHD0400009897
|
BTA4
|
MY,MC,R,M/C,H,P,M
|
SW/GR
|
0.57
|
BovineHD0500015340
|
BTA5
|
MC,R,M/C,H
|
SW/VC
|
0.69
|
BovineHD2900011510
|
BTA29
|
ALL TRAITS
|
SW/HR
|
0.53
|
BovineHD0300002536
|
BTA3
|
MC,P
|
SW/KG
|
0.80
|
BovineHD2400003366
|
BTA24
|
R,M/C,H,M
|
SW/OG
|
0.64
|
BovineHD0800001684
|
BTA8
|
R,M/C,P
|
TP/GR
|
0.59
|
BovineHD0300009588
|
BTA3
|
ALL TRAITS EXCEPT M
|
TP/VC
|
0.73
|
BovineHD1900008163
|
BTA19
|
ALL TRAITS EXCEPT MY
|
TP/HR
|
0.50
|
BovineHD2800011107
|
BTA28
|
M/C
|
TP/KG
|
0.81
|
BovineHD1300018382
|
BTA13
|
MC,R,M
|
TP/OG
|
0.62
|
BovineHD0200032623
|
BTA2
|
ALL TRAITS
|
GR/VC
|
0.81
|
BovineHD1900018844
|
BTA19
|
ALL EXCEPT MY
|
GR/HR
|
0.60
|
BovineHD1000011926
|
BTA10
|
ALL EXCEPT MY,M
|
GR/KG
|
0.81
|
BovineHD0500019637
|
BTA5
|
ALL TRAITS
|
GR/OG
|
0.65
|
BovineHD0500012675
|
BTA5
|
ALL TRAITS
|
VC/HR
|
0.80
|
BovineHD1900008163
|
BTA19
|
ALL EXCEPT MY
|
VC/KG
|
0.76
|
BovineHD0200022436
|
BTA2
|
R,M/C,H,P
|
VC/OG
|
0.61
|
BovineHD1200026864
|
BTA12
|
ALL EXCEPT MY
|
HR/KG
|
0.80
|
BovineHD1300013278
|
BTA13
|
ALL TRAITS
|
HR/OG
|
0.60
|
BovineHD1700005865
|
BTA17
|
ALL EXCEPT MY,H
|
KG/OG
|
0.78
|
BovineHD2100009715
|
BTA21
|
R,M/C,H,P
|
*MY-Milk Yield, MC-Milk Composition, R-Reproduction, M/C-Meat & Carcass |
H-Health, P-Production, M-Morphology |
**SW-Sahiwal, TP-Tharparkar, GR-Gir, VC-Vechur, OG-Ongole, KG-Kangayam, HR-Hariana |
The annotation of 5 top most genomic regions among these cattle to Bovine QTL database revealed that these regions harbored QTL affecting milk yield, milk composition, reproduction, production, Health status and morphological traits (Supplementary Table S2). The top divergent regions among most divergent breed pairs may be putative selection signature for diversifying traits between the breeds. For example, the marker BovineHD0500014902 lies in the FAM19A2 gene and is the most differentiating region between Kangayam and Sahiwal. This gene is thought to produce brain-specific chemokines or neurokines, which function as immune and nervous cell regulators (https://www.genecards.org/cgi-bin/carddisp.pl?gene=TAFA2), and thus may affect Kangayam and Sahiwal's health in different ways. Kangayam is a hardy draft breed compared to Sahiwal. There were significant differences in number of QTLs affecting milk yield and production traits when Kangayam paired with any of the dairy breed. For example, between Hariana and Kangayam, the marker BovineHD0500012581 lies in RAB31 gene which is involved in metabolic pathway (https://www.genecards.org/cgi-bin/carddisp.pl?gene=RAB31) and has also been annotated with QTLs affecting milk performance traits (Cattle QTL data base). Hence, may be responsible for variations in these traits.
The other selection signatures identified between Kangayam and dairy breeds included several genes like Bestrophin 3 (BEST3) [https://www.genecards.org/cgi-bin/carddisp.pl?gene=BEST3&keywords=BEST3], Diacylglycerol Kinase Alpha (DGKA) [https://www.genecards.org/cgi-bin/carddisp.pl?gene=DGKA&keywords=DGKA], Adenosylhomocysteinase (AHCY) [https://www.genecards.org/cgi-bin/carddisp.pl?gene=AHCY&keywords=AHCY], Phosphatidylinositol Glycan Anchor Biosynthesis Class U (PIGU) [https://www.genecards.org/cgi-bin/carddisp.pl?gene=PIGU&keywords=PIGU] and Phosphofructokinase, Platelet (PFKP) [https://www.genecards.org/cgi-bin/carddisp.pl?gene=PFKP&keywords=Phosphofructokinase] which are involved in transportation of glucose and sugars, signaling pathways, cellular processes, cell division and glycolysis regulation, respectively. Moreover, these genomic regions also harbour QTLs affecting milk performance traits (Cattle QTL data base). The signatures were also identified even between the dairy breeds. For example, the genomic region, having BovineHD1300006213 locus, revealing high differentiation between Tharparkar and Gir (FST=0.585) harbours QTLs affecting all the traits under study except milk yield (Cattle QTL data). There were significant differences in number of QTLs affecting production (body weight, growth rate etc.) and morphological traits (height etc.) of Vechur when paired with large sized cattle (Supplementary Table S3). The genomic region having BovineHD0500013882 locus starting from 48008400 to 48069099 on chromosome 5 is having HMGA2 gene which is responsible for explaining the variation in stature both in cattle42 and human beings43.
Implications Of Genomic Analyses For Breed Management And Conservation
The genomic analyses detailed herein revealed that all the high yielders dairy breeds viz. Sahiwal, Tharparkar and Gir shared common ancestry but are admixed due to gene flow among them through crossing, migration and /or grading up with other breeds. Hariana, traditionally a dual purpose breed but was more closely associated with dairy breeds. Ongole, Kangayam and Vechure were quite distinct from rest of the breeds under study. The genomic analyses represented the unique gene pool of these cattle genetic resources, befitting their breed’s status. The Hariana cattle was under intensive selection for high milk yield in the past and was being used as improver breed for many other breeds/populations for augmenting milk production in the country, and same was also revealed here through genomic characterization.
In conclusion, BovineHD BeadChip genotyping of Indian cattle breeds is promising, as it can not only structure these cattle but also detect distinct signatures of selection. As a result, it could be used for a wider range of studies, such as genome wide association studies and genomic selection involving larger populations of the breeds. Multi-breed genomic selection may be feasible in dairy breeds due to their shared genome. This study reveals a trend toward shrinking effective population sizes in native Indian cattle breeds, indicating that a long-term breeding strategy is needed to prevent further reductions in Ne, as well as genetic improvement and potential conservation. In future, whole genome sequencing information on these breeds may be useful for pinpointing the genomic regions linked to polygenic productivity, health, fertility, and behavioural traits that evolved under Indian ecological and farming systems.