Diversity assessment of a lesser known buffalo population from Central India and its comparative evaluation reveals presence of sufficient genetic variation and absence of selection

Abstract

Planned breeding and conservation strategies for a lesser-known population require an assessment of complete genetic diversity and population structure analysis in addition to its morphometric characteristics. In the present study, a comparative analysis of the genetic structure of a rare buffalo population, namely Chhattisgarhi, was extensively studied using a panel of FAO-recommended microsatellite markers along with well-established breeds namely Murrah, Nili-Ravi, Gojri, Kalahandi, and Nagpuri. Mode shift analysis indicated the absence of genetic bottleneck in the recent past. Assessment of genetic diversity indices across all loci indicated the presence of sufficient genetic variation within and between populations. Analysis of molecular variance between the six different buffalo populations attributed 19.05% of the variations to between-population differentiation. Cluster analyses using DAPC and Bayesian approach along with the phylogenetic tree based on UPGMA grouped six populations into three groups. The Chhattisgarhi population was revealed to be genetically closer to Nagpuri and Kalahandi populations. The study reveals the presence of sufficient genetic diversity within the Chhattisgarhi population and indicates the absence of a systematic selection program. We suggest improvement and conservation programs should be planned for this breed in the near future through short-term selection.

Keywords:

• Chhattisgarhi buffalo
• genetic diversity
• microsatellite markers

Introduction

In the recent past, a large population of interbred non-descript buffaloes have evolved¹ causing a negative impact on the diversity and uniqueness of comparatively lesser-known local breeds of buffalo like Chhattisgarhi. As per the² 79.7% of buffaloes are either non-descript or grades of pure breeds. A genetic marker-based examination of a population structure and its comparison with breeds from the surrounding area and well-known milch buffalo breeds can be used to genetically characterize a population. Being neutral, codominant, ubiquitous, and highly polymorphic, microsatellite markers have been most commonly used for the study of genetic diversity in different domestic animal species.³ In the present study, comparative assessment of the genetic structure of Chhattisgarhi buffalo population with other well-established breeds of buffaloes was done using 15 microsatellite markers. Chhattisgarhi buffalo has its breeding tract in the northern hilly and central plain areas of Chhattisgarhi (Figure 1) state and has been defined well based on its morphological characteristics and is reared under unique system of management ‘Bathan’.⁴ Five different buffalo breeds namely Murrah, Nili-Ravi, Gojri, Nagpuri, and Kalahandi were used for comparative analysis. Such a comparative genetic study identifies the hidden production potential of the uniquely adapted breeds to various geographical regions of India, and their genetic characterization will aid in formulating breed-specific conservation and genetic improvement programs.

Figure 1. Map highlighting the breeding tract of Chhattisgarhi buffalo.

Materials and methods

Sampling and DNA extraction

With the consent of herd owner, 5-10 ml of blood samples of 48 healthy buffaloes representing the original indigenous true-to-type phenotype of Chhattisgarhi buffalo population were selected for this study from the native tract following the guidelines of measurement of the domestic animal diversity program.⁵ Whole genomic DNA was isolated using the standard phenol-chloroform method described by Sambrook and Russel.⁶ The primers for 24 microsatellite markers as suggested by ISAG-FAO were used. The genotyping was carried out by following the protocol of Mishra et al.⁷ and Singh et al.⁸ with ABI PRISM 3100 automatic sequencer (Applied Biosystems, Foster City, CA, USA). To maximize the throughput, products were amplified by different primers with different dyes were pooled for one capillary injection. GeneScan-500 LIZ™ Size Standard (Applied Biosystems) was used as the internal standard for fragment sizing. Internal lane size standard was run with every sample for accurate sizing. Gene Scan software was used to calculate the molecular weight (size) of different alleles/DNA fragments.⁹ For comparative analysis 15 microsatellite markers (CSSM45, CSSM47, CSSM66, HEL13, ILSTS19, ILSTS25, ILSTS26, ILSTS28, ILSTS29, ILSTS30, ILSTS33, ILSTS56, ILSTS58, ILSTS60, ILSTS61) common in all six populations were taken (data for Kalahandi, Murrah and Nili-Ravi populations was taken from Animal Genomics lab, NBAGR, Karnal).

Statistical analysis

BOTTLENECK v1.2.01¹⁰ was used to test for the departure from the mutation drift equilibrium and the presence of any bottleneck based on heterozygosity excess and heterozygosity deficit. The variance and probability were set at 12.00 and 95%, respectively for the estimation of Bottleneck with 5,000 simulations. Data analysis was done on¹¹ platform with the help of packages such as ‘adegenet’ v2.1.7¹² ‘poppr’ v2.9.3¹³ and ‘hierfstat’ v0.5-11¹⁴ ‘polysat’ v1.7¹⁵ and ‘PopGenReport’.¹⁶ Computation of genetic diversity indices, allele frequencies, and AMOVA were performed through the mentioned R packages. The effective number of alleles (N_E) was calculated by using the software GenoDive©2005-2022.¹⁷ UPGMA¹⁸ approach was used to construct a phylogenetic tree using the package ‘ape’ v5.6-2.¹⁹ DAPC (Discriminant analysis of Principle components analysis) was performed with the package ‘ada4’²⁰ and ‘adegenet’ v2.1.7.¹² Population assignment was done by the ‘ParallelStructure’ package, which is based on STRUCTURE v2.3.4.²¹ The most likely number of subpopulations was found by the Evanno ΔK method²² using the package ‘POPHELPER’ v3.0.4.²³

Results

The status of any bottleneck in the population was inferred from the sign, standardized, and Wilcoxon tests. There were 14 loci with significant heterozygosity excess for all three models (IAM, TPM, and SMM) (Table 1). Probability values for the null hypothesis for bottleneck were calculated through Mode shift analysis that followed L-shaped curve when microsatellite alleles were classified into ten frequency classes (Figure 2).

Figure 2. L-shaped mode shift curve showing absence of recent genetic bottleneck in Chhattisgarhi population.

Table 1. Mutation-drift equilibrium, heterozygosity excess under different mutation models in Chhattisgarhi population.

Test	Model
		IAM	TPM	SMM
Sign test (Number of loci with heterozygosity excess)	Exp	14.06	14.14	14.14
	p-value	0.00936	0.39146	0.00003
Standardized differences test	T2 value	1.882	−2.230	−9.346
	p-value	0.02991	0.01288	0.00000
Wilcoxon rank test (one tailed for heterozygosity excess	p-value	0.00752	0.61585	0.99990

IAM: infinite allele model; TPM: two phase model; SMM: stepwise mutation model.

Locus-wise details of overall populations are given (Table 2). No. of alleles across all the loci observed in the study were 270. A total of 112, 128, 160, 103, 91, and 92 alleles were found in Chhattisgarhi, Gojri, Kalahandi, Murrah, Nagpuri, and Nili-Ravi populations, respectively. Population-wise values for comparison among the populations are given (Table 3). Analysis of molecular variance (AMOVA) revealed 19.05% of total variations were present between studied populations, while 74.14% of variations in the population were within individuals of the population. A 6.80% of variation was attributed to between individuals within the population. The pair-wise F_ST values were used to illustrate the genetic distance between populations, which significantly differentiated all six buffalo populations, as shown in the F_ST heat map (Figure 3). Chhattisgarhi buffaloes were found to have a low genetic distance from Nagpuri and Kalahandi while a higher genetic distance from Murrah and Nili-Ravi. UPGMA tree plot results showed two main clusters among all studied populations (Figure 4). Discriminant analysis of principal components revealed that the first discriminant function on principal components explained 44.8% of the genetic variance. Populations in the data were grouped into three clusters as represented in the scatter plot (Figure 5). Population assignment was accomplished in STRUCTURE for K = 2 to 7 and results are presented in the form of a bar plot (Figure 6). Evanno ΔK method suggested the optimum value for K was 3 (Figure 7). The proportion of the membership of each pre-defined population in each of the 3 clusters at a K value of 3 is given (Table 4).

Figure 3. Heat-map showing pairwise F_ST among six buffalo populations.

Figure 4. UPGMA tree inferring evolutionary distances of the six buffalo populations.

Figure 5. Scatter plot of individuals from buffalo populations clustered through discriminant analysis of principal components (DAPC).

Figure 6. Cluster of individuals from different populations through bayesian genetic approach of STRUCTURE. Bars depict the individuals in a population and the color corresponds to the inheritance from different breeds. The 6 population studied are represented as CGH – Chhattisgarhi, GOJ – Gojri, KLH – Kalahandi, MUR – Murrah, NAG – Nagpuri, and NIR – Nili-Ravi.

Figure 7. Evano’s plot showing the most likely number of genetic clusters present in the population.

Table 2. Locus-wise genetic indices based on overall microsatellite data of six buffalo populations.

Locus	NA	NE	Ho	HE	FST	FIS	PIC
CSSM45	14	3.675	0.637	0.738	0.147	0.137	0.867
CSSM47	36	4.715	0.679	0.799	0.143	0.151	0.932
CSSM66	25	3.688	0.581	0.739	0.177	0.214	0.893
HEL13	13	4.157	0.786	0.768	0.1	−0.024	0.838
ILSTS19	12	1.286	0.229	0.225	0.601	−0.02	0.481
ILSTS25	16	3.059	0.672	0.681	0.133	0.013	0.763
ILSTS26	13	2.502	0.634	0.607	0.238	−0.044	0.779
ILSTS28	17	3.995	0.751	0.759	0.006	0.01	0.740
ILSTS29	19	2.163	0.545	0.544	0.341	−0.002	0.821
ILSTS30	18	3.072	0.669	0.683	0.197	0.02	0.846
ILSTS33	12	3.048	0.547	0.682	0.117	0.197	0.790
ILSTS56	21	2.410	0.524	0.593	0.221	0.115	0.772
ILSTS58	16	5.551	0.832	0.83	0.043	−0.003	0.852
ILSTS60	19	2.465	0.468	0.604	0.169	0.225	0.756
ILSTS61	19	3.818	0.611	0.749	0.095	0.183	0.826
Mean	18	3.807	0.611	0.667	0.182	0.078	0.797

Table 3. Population-wise genetic indices of six buffalo populations.

Population	N	NA	MNA	AR	Ho	HE	Fis
CGH	48	112	7.467	6.928	0.607	0.668	0.088
GOJ	48	128	8.533	7.93	0.783	0.735	−0.074
KLH	48	160	10.667	9.876	0.606	0.753	0.188
MUR	40	103	6.867	6.733	0.524	0.608	0.137
NAG	48	91	6.067	5.681	0.598	0.63	0.045
NLR	40	92	6.133	5.976	0.549	0.608	0.093

Table 4. Results of STRUCTURE showing inferred cluster and ancestry of buffalo populations at K = 3.

Given population	Given Inferred Clusters			Number of individuals
	1	2	3
CGH	0.005	0.005	0.990	48
GOJ	0.990	0.003	0.007	48
KLH	0.006	0.012	0.982	48
MUR	0.010	0.980	0.010	40
NAG	0.003	0.005	0.993	48
NLR	0.002	0.994	0.004	40

Discussion

For the Standardized differences test, the p-value for the null hypothesis for mode shift equilibrium was significant for both TPM and SMM models. For the Wilcoxon test, which is a more robust and commonly used test for bottleneck analysis than the other two tests,²⁴ the p-value for the null hypothesis for mode shift equilibrium in both TPM and SMM models is > 0.05 thereby indicating no genetic bottleneck in recent past. Graphical representation of Allelic classes and proportion of alleles²⁵ was also observed to be a normal L-shaped distribution suggesting no evidence of a severe reduction in effective population size in the immediate past, thus no bottleneck.

The mean H_E of 0.667 across all loci and populations indicated sufficient polymorphism to measure genetic variation.²⁶ The mean N_A per locus was 7.467 for Chhattisgarhi; these were highest in Kalahandi (10.67) and lowest in Nagpuri buffaloes (6.06). The mean no. of alleles per locus in the present study is higher than reported by Mishra et al.²⁷ in Bunny and Murrah populations and by Vijh et al.²⁸ in 12 buffalo populations, however, El-Kholy et al.²⁹ found a higher Mean no. of alleles per locus among three buffalo breeds in Egypt. The mean expected heterozygosity (0.71-0.78) in the study was in the range of 0.608 to 0.753. Kumar et al.³⁰ reported higher mean expected heterozygosity (0.71-0.78) in a study involving eight populations, whereas²⁷ reported lower values (0.57-0.659). In the present study, buffalo populations from different regions depict a significant differentiation, attributable to their genetic distinctiveness. An overall moderate amount of inbreeding in the population was indicated by heterozygote deficiency was evident from the positive mean F_IS value. This could be possible due to the mating of close relatives in the herds owned by livestock owners or due to genetic hitchhiking.³¹ 19.05% of variations between populations from results obtained from AMOVA suggest the distinctness of all six breeds. Between populations, variations were almost double that²⁷ reported among five buffalo populations. As shown in the pair-wise F_ST heat map Chhattisgarhi buffalo was least differentiated from Nagpuri and Kalahandi, moderately from Gojri and the highest differentiation was found from Nili-Ravi and Murrah based on F_ST. Mishra et al.²⁶ and Kathiwavan et al.³² reported F_ST values of 0.145 and 0.083 among the buffalo populations studied, which were lesser than the overall F_ST value of 0.173 in the present study. These results were also in compliance with the results from the factorial correspondence analysis based on molecular data and phylogenetic tree obtained from the dissimilarity matrix. According to UPGMA tree diagram (Figure 4) with two central clusters I and II, Chhattisgarhi and Nagpuri were grouped under one group (Iaa) and Kalahandi as a separate subgroup (Iab) in the same sub-cluster, while Gojri was placed under the cluster Ib. Murrah and Nili-Ravi were completely clustered differently as second central cluster (II). In the scatter plot of DAPC, Chhattisgarhi buffaloes are clustered with Nagpuri and Kalahandi, whereas Murrah and Nili-Ravi are clustered together. Meanwhile, Gojri was plotted on the opposite side of axis 2. STRUCTURE analysis for K = 2 shows there was no differentiation between Chhattisgarhi, Gojri, Nagpuri, and Kalahandi populations, whereas Murrah and Nili-Ravi populations were grouped under one cluster. For K = 3, Chhattisgarhi, Nagpuri, and Kalahandi buffaloes were grouped under one cluster, Gojri buffaloes as a separate cluster while Murrah and Nili-Ravi were grouped under the third cluster. Similar results for Murrah, Nili-Ravi and Gojri populations were reported by Kathiwavan et al.³³ This indicated that studied populations reared in other geographic regions have differentiated and possess unique allelic combinations than those reared in similar geographical areas. Presently the breeding areas of all these similarly grouped populations are overlapping.

Conclusion

Microsatellite loci used in the study were found sufficiently polymorphic, indicating that they were suitable for genetic diversity analysis in buffalo populations. Chhattisgarhi buffalo population is genetically closer to the Nagpuri and Kalahandi buffalo breeds but well differentiated from the Murrah and Nili-Ravi breeds of buffaloes. Hence, based on the findings of the present study it is recommended to State or Central Govt. implementing agencies to formulate an improvement and conservation program, which in long run will establish this breed and shall uplift the sustenance income of the local buffalo-rearing farmers from Central India.

Acknowledgments

The authors thank the Directors of ICAR-NDRI and ICAR-NBAGR for their support to carry out the study. Technical support received from Mr. Subhash Chander, T-5, NBAGR, is gratefully appreciated.

Disclosure statement

All authors declare that they have no conflict of interest.

Additional information

Funding

The author(s) reported there is no funding associated with the work featured in this article.