Detection of Selection Signatures in some of the Water Buffaloes across the World

doi:10.21203/rs.3.rs-4516365/v1

Download PDF

Research Article

Detection of Selection Signatures in some of the Water Buffaloes across the World

https://doi.org/10.21203/rs.3.rs-4516365/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

In order to identify the selection signatures of the water buffalos across the world, the genomic information of 165 buffalos which belonged to 15 genetic groups of buffaloes was used. The genomic information was obtained from Dryad (doi:10.5061/dryad.h0cc7). The quality control and data filtration were performed using PLINK1.9 software. The genetic clustering and the population structure was examined using the GenABEL and Admixture1.23 software's, respectively. The results of principal component analysis showed that the examined populations could be classified into 4 separate categories. The results of population structure analysis confirmed the results of principal components analysis. The signatures of selection were searched with the help of iHS statistics using the ReHH software. Moreover, the unbiased F_ST (θ) estimator was calculated using the Plink1.9 software. The 25 and 24 genomic regions, which passed the unbiased F_ST and iHS statistics thresholds, were identified as selection cues, respectively. Selected regions were aligned on the bovine genome and 411 genes related to selected regions were identified. Of all the identified genes, 53 genes related to olfactory receptors (OR), 51 genes somehow involved in cell membrane structure and animal immunity against pathogens including initiate and regulate the immune response. The identified QTLs related to detected regions, were associated with milk production, milk somatic cells, fertility, ion disease, calving and growth. There is an acceptable consistent between the milk and fat production genes and the related identified QTLs.

Buffalo

Genome

Signatures of selection

In general, obtaining information on the origin of the breed and the history of the evolution of populations is important to predict the gene structure of each breed and to describe different traits such as disease resistance, stress tolerance, and adaptation to different environments in the future [1]. The results of the studies of domestic livestock have shown that the genetic diversity of domestic animals has decreased during the domestication process due to various genetic factors such as genetic drift, mutation, natural and artificial selection of breeds, and the use of a small number of breeds for more economic production. The decrease in the diversity of these animals has caused global concern [2]. In this regard, using genome information in order to identify the structure and the evolutionary history of populations, to develop strategies to conserve the diversity of breeds and to livestock breeding program is very practical and useful. The identification of the signatures of selection is regarded to be one of the important and practical aspects of population genomic analysis. The signatures of selection refer to the patterns on the genome which are formed on specific regions of the genome by the selective forces [3]. They indicate the decrease in the spatial genetic diversity, the deviation from the site frequency spectrum, the increase in linkage disequilibrium and the development of haplotype [4]. Different tools and methods have been developed and used successfully in order to identify the signatures of selection at the genome level of various populations. In general, these methods can be classified into 5 distinct categories including statistical tests based on population differentiation [5], statistical tests based on linkage disequilibrium and haplotype length [6], statistical tests based on site frequency spectrum [7], statistical tests based on the decrease in spatial genetic diversity [8], and the statistical tests based on performance-altering mutations [9]. Moreover, there are a number of tests which constitute a combination of methods for identifying the signatures of selection. The results of the relevant studies have highlighted the relative sensitivity of combined methods for identifying the signatures of selection [10]. Simianer et al. (2014) examined the efficiency of the most important and the most common methods of identifying the signatures of selection [3]. The domestic animals are regarded to be more appropriate cases regarding the identification of the signatures of selection in comparison with humans and are used to identify the genes that control different phenotypes. This issue stems from the fact that, in recent years, these animals have been exposed to intense artificial and natural selection [11]. To date, studies have been carried out for identify genomic regions in several livestock including cattle [12–17], goats and sheep [18–20] and pig [21, 22]. Similar studies have been performed on other livestock species genome including buffalo [23–25], horses [26, 27], camel [28, 29] and chiken [8, 30, 31]. In most studies, detected genes are related to adaptation and performance traits. The present study intends to use the genomic information of water buffaloes in buffalo breeding countries in order to identify the regions of the buffalo genome that have been affected by natural or artificial selection forces for many years. These genomic regions may be the genomic regions which distinguish the above-mentioned breeds. The results of the study can be useful for analyzing the differences between the examined buffalo breeds and implementing breeding programs.

Genomic Data

The present study used the genomic data of 165 water buffaloes which belonged to 15 breed populations and were from 11 major buffalo breeding countries. The genomic information of water buffaloes across the world was obtained from Dryad information repository (doi: 10.5061 / dryad.h0cc7). The used samples accounted for a large part of the geographical distribution of the water buffalos across the world. The sequenced samples were obtained from 11 countries, including Italy, Mozambique, Romania, India, Bulgaria, Brazil, Turkey, Egypt, Iran, Pakistan, and Colombia [32]. The Mediterranean and the Murrah breeds are among the most prominent breeds of buffalo in the world. Most of the breeding programs have focused on these breeds. At the present time, these breeds are being exported to the other countries. The samples, which were used in this study, were provided by the members of the International Water Buffalo Consortium. All of the samples were sequenced in the Affymetrix Laboratory (Santa Clara CA USA) using buffalo-specific Axiom R Buffalo Genotyping Array (http://www.affymetrix.com).

Quality Control and Data Screening

Quality control and filtering was performed using the PLINK software [33]. Genomic data of 56,845 SNPs from 165 animals were qualitatively controlled as follow thresholds. First, the individuals and SNP with more than 5% missing genotype, were excluded. Then, the data were screened based on Minor Allele Frequency (MAF) less than 1%, and Hardy-Weinberg Equilibrium (HWE).

Principal Component Analysis and Population Structure

At this stage, the Principal Component Analysis (PCA) and population structure analysis were performed to examine the genetic diversity and to determine the genetic groups. PCA is used as the first analysis of data investigation and data description in most population genetic analyses [34]. The PCA analysis was performed using GenABEL package [35] in R software, and the population structure analysis was carried out using Admixture 1.23 software [36].

Detection of Selection Signatures

The signatures of selection were searched unbiased fixation index (F_ST) using PLINK1.9 software, and iHS statistics with ReHH software [37] after being imported and phased using Beagle software [38].

Biological Pathway Analysis

The genomic regions which were selected using the signatures of selection identification method were aligned on cow genome (ARS-UCD1.2) on the online Ensemble Biomart Tool database (http://asia.ensembl.org/biomart/martview). Next, all of the genes, which were associated with these genomic regions, were identified. At the present time, this database contains information on 27607 genes. At this stage, first, the selection sites of the signatures of selection search methods were annotated with the location of genes which were listed in the buffalo genome in the Ensemble Biomart Tool. Second, the list of the genes, which were related to the selected regions, was obtained. Identifying biological pathways and gene networks was done by DAVID 2021 (The Database for Annotation, Visualization and Integrated Discovery) [39]. Gene interaction networks and functional enrichment analysis was performed using STRING 11.5 [40].

Quality control and data screening

Among all of the available individuals and markers, 10 individuals were excluded due to missing genotype. Moreover, 212 SNP markers were excluded because of missing genotype of more than 5%. Finally, 4416 SNP markers were excluded owing to MAF of more than 1%. In addition, 576 SNPs were excluded because they were not in Hardy-Weinberg equilibrium. The threshold level of Hardy-Weinberg equilibrium was 8.8 ×10 ^− 7. The genotypic rate in 51655 markers of the remaining 155 individuals was equal to 0.9973. Moreover, 18 markers were excluded from among the markers which had acceptable quality and were not examined in the further analysis due to their unknown genomic location. Moreover, 1595 markers were excluded since they were located on the sex chromosome. Lastly, 50042 SNP markers were identified for population structure analysis using F_ST-based selection criteria and iHS-based analysis. The distribution of markers that passed the quality control steps is shown in the Fig. 1 in appendix based on chromosomes (Supplementary Fig. 1).

Principal Component Analysis (PCA) and population structure

The results of the principal component analysis and diversity studies showed that the examined populations could be classified into four separate categories (Fig. 1 and Supplementary Fig. 2).

[Figure 1 Position]

Figure 1. PCA analysis of 15 the world river buffalo populations using genomic information

Moreover, the results of the population structure analysis, which was performed using Admitxure software, confirmed the results of the principal component analysis. Similarly, the results of this section underlined the existence of three genetically differentiated groups with higher purity and a mixed group of the existing genetic bases (Fig. 2). The first group included the Italian, Mozambican, and Mediterranean breeds. The second group involved a combination of the other breeds and constituted a separate group. The third group was comprised of Pakistani, Indian, Brazilian-Murrah, Bulgarian-Murrah, and Colombian breeds. Finally, the fourth group included Iranian, Egyptian and Anatolian (Turkish) breeds. The results of the population structure analysis, which was performed using Admitxure software, confirmed the results of the principal component analysis.

[Figure 2 Position]

Figure 2. Graph of population structure and admixture of 15 the world river buffalo populations. Each vertical column belongs to one individual. Each column was represented by one or combination of several colors. The every color in the structure of individual's genome indicate a specific genetic origin

Detection of selection signatures

In the F_ST analysis, which was performed to identify the signatures of selection, 24 regions of the genome were identified as signatures of selection. The numerical theta value of these regions was higher than 0.22 and they included 0.1% of the examined SNP markers. The identified regions were located on chromosomes 2, 4, 6 (3 regions), 8 (4 regions), 9, 10, 11 (3 regions), 13, 17, 18 (3 regions), 23 (2 regions), and 26 and 27 (2 regions) respectively (Fig. 3). The analysis showed that, only 17 regions of the above-mentioned regions contained the gene.

[Figure 3 Position]

Figure 3. The Manhattan plot of the 15 river buffalo populations overall unbiased F_ST based on averaged values with 300Kbp window length. The red line indicates the threshold for selection signs (on autosomal BTA). The threshold in this study was determined based on the experimental F_ST distribution.

The iHS statistic is a powerful tool for examining LD destruction around the selected region on genomic by assessing haploid characteristics within a population [6]. The regions with high frequency and high iHS, which are near the selected allele, are the targets of positive selections. In the case of cattle, which are selected both naturally and artificially, these regions can be the candidates for the major effect genes [12]. In the examination of the signatures of selection using iHS method, 25 genomic regions, which exceeded the threshold level of this statistic (The numerical value of this statistic, which was calculated by logarithmizing the determined numerical values, was equal to 4), were identified as signatures of selection. The identified regions were located on chromosomes 1, 2, 3, 4, 5 (2 regions), 6 (2 regions), 8, 9, 10 (2 regions), 12 (2 regions), 15 (4 regions), 16, 18 (2 regions), and 20, 27 and 29 (2 regions) respectively (Fig. 4). Moreover, 184 genes were identified in these 15 selected regions. The regions, which contained genes, were located on chromosomes 1, 2, 3, 4, 5, 6, 9, 10 (2 regions), and 12, 15, 18, 20 and 29 (2 regions) respectively. The selected region on chromosome 5 was one of the important identified regions in this study, because it contained many genes and was vital. This region and the regions of chromosomes 10 (position at 25 Mb), 12 (position at 60 Mb) and 15 (position at 60 Mb), and the related genes were identified as the signatures of selection in the comparison between Azeri and Khuzestani buffalo breeds [41].

[Figure 4 Position]

Figure 4. The Manhattan plot of iHS statistics of the 15 river buffalo populations. The red line indicates the threshold for selection signs (on autosomal BTAs).

QTLs related to the selected genomic regions

The 227 genes, which were related to 17 of selected regions from averaged unbiased fixation index (F_ST), were identified. Based on the analysis, 30 of these genes were related to the olfactory receptors and about 50 genes were associated with cell membrane and immune response in animals. Moreover, 28 genes of the genes, which were related to the olfactory receptors, were located in the range of 29.16 to 29.93 Mb of the region of chromosome 23 and 2 genes were located on chromosome 4. These genes included OR11A1, OR11A12, OR12D3, OR12D20, OR12D20, OR10AL40 and OR10AL41 among others. In the study by Mokhber et al. (2018), the 1 Mb selected region on chromosome 23 was identified as a signature of selection [23]. This region is rich in olfactory-related genes. Fifty-one of the identified genes were related to the cell membrane, influenced animal immunity against pathogens, and were involved in the initiation and regulation of the immune responses in some way. These genes were located on chromosomes 4, 6, 17 and 23. Most of the genes, which were identified in this region, were related to chromosomes 17 and 23 in terms of density.

Thirty-on genes of the 184 identified genes from iHS statistic, included 4 cases of lncRNA, 5 cases of miRNA, 3 cases of misc_RNA, 8 cases of snRNA, 7 cases of snoRNA, and 4 cases of non-protein-coding pseudo-genes. The remaining 153 genes were related to the genes encoding proteins. Twenty-three of these codes were related to the olfactory receptors (ORs) which are commonly located on chromosomes 10 and 15. They included OR6S1, OR2D2, OR6A2, OR52B2, OR56A1, OR56A5 and OR52N1 among others. The olfactory receptor genes have been reported as signatures of selection in humans [42], domestic animals such as dogs [43], pigs [44], cows [4] and buffalos [23]. The olfactory gene family constitutes a group of the genes which are involved in the evolution and domestication.

[Table 1 position]

Table 1

Complete list of genomic regions and genes harboring significant SNPs identified by unbiased FST and his methods
Chr	Start	End
1	93.39	93.40	iHS	ENSBTAG00000046279, ENSBTAG00000012401
2	10.07	11.07	FST	ENSBTAG00000049344
2	76.11	76.92	iHS	CNTNAP5
3	28.71	29.69	iHS	DENND2C, BCAS2, TRIM33, SYT6, AP4B1, BCL2L15, PTPN22, RSBN1, PHTF1
4	94.06	95.23	iHS	CPA5, CPA1, ENSBTAG00000054642, ENSBTAG00000052730, ENSBTAG00000048995, MEST, bta-mir-335, ENSBTAG00000051746, MKLN1
4	105.89	106.88	FST	ENSBTAG00000050190, ENSBTAG00000053376, TRBV24-1, ENSBTAG00000053701, ENSBTAG00000053785, ENSBTAG00000054917, ENSBTAG00000051147, TRBV29-1, PRSS2, LOC509513, LOC100300510, EPHB6, OR(ENSBTAG00000048380), OR6V1, PIP ENSBTAG00000050494, TAS2R39, TAS2R40, GSTK1, TMEM139, CASP2, CLCN1, ZYX, TAS2R60, ENSBTAG00000031162
5	55.63	57.78	iHS	CTDSP2, bta-mir-26a-2, AVIL, METTL1, CYP27B1, CDK4, AGAP2, ENSBTAG00000050652, B4GALNT1, ENSBTAG00000051593, DCTN2, DDIT3, ARHGAP9, R3HDM2, STAC3, NDUFA4L2, SNORA62, STAT6, NEMP1, MYO1A, TAC3, ZBTB39, ENSBTAG00000049081, ENSBTAG00000050051, RDH16, SDR9C7, ENSBTAG00000055155, ENSBTAG00000042169, ENSBTAG00000043104, PRIM1, NACA, PTGES3, ATP5F1B, ENSBTAG00000043048, bta-mir-677, bta-mir-677, BAZ2A, GLS2, MIP, TIMELESS, APOF, APON, STAT2, PAN2, CNPY2, bta-mir-12054, CS, ANKRD52, ENSBTAG00000052361, RNF41, SMARCC2, PMEL, PYM1, MMP19, DNAJC14, SARNP, ENSBTAG00000051365, CD63, ITGA7, OR10P25, ENSBTAG00000054814, ENSBTAG00000052662
6	1.34	1.65	FST	ENSBTAG00000051456, MARCHF1
6	83.29	83.35	iHS	STAP1
6	88.51	90.6.8	FST	AFP, AFM, ENSBTAG00000049436, 7SK, ENSBTAG00000051518, ENSBTAG00000052772, CXCL8, ENSBTAG00000027534, CXCL5, PF4, CXCL2 ENSBTAG00000051891, CXCL3, GRO1, MTHFD2L, EPGN, EREG, AREG, PARM1, THAP6, ENSBTAG00000050403, ENSBTAG00000046788, USO1, U1
8	7.23	7.92	FST	ENSBTAG00000053818, ENSBTAG00000025954, ENSBTAG00000008678, ENSBTAG00000039873DEFB134, DEFB136, CTSB, FAM167A
8	19.63	19.95	iHS	IZUMO3, U6
8	26.11	26.12	FST	ENSBTAG00000052821
8	55.57	55.72	FST	TLE4
8	77.95	79.25	FST	NTRK2, ENSBTAG00000050767, U6, NAA35
9	48.47	48.48	iHS	ENSBTAG00000033109
9	81.65	83.15	FST	UTRN, ENSBTAG00000054706
10	25.53	26.56	iHS	OR4E1, OR10G1, OR10G1C, OR10G3, SALL2, METTL3, RAB2B, CHD8, ENSBTAG00000042579, ENSBTAG00000043196, SUPT16H, HNRNPC, OR5AU1, ZNF219, RNASE13, NDRG2, ENSBTAG00000049034, RNASE1, ENSBTAG00000053798, BRB, OR6S1, RNASE12, RNASE11, ENSBTAG00000003006, ENSBTAG00000052633
10	37.34	38.34	iHS	GANC, CAPN3, SNAP23, HAUS2, ENSBTAG00000048958
10	86.49	87.37	FST	FOS, JDP2, bta-mir-10162, BATF, FLVCR2, ENSBTAG00000011985, TTLL5, IFT43
11	67.01	67.44	FST	PROKR1, ARHGAP25, GKN3P, GKN1, ANTXR1, 7SK
11	70.15	71.32	FST	ALK, PCARE, ENSBTAG00000050766, TRMT61B, ENSBTAG00000048347
12	44.56	44.58	iHS	OR6A2
12	62.29	62.30	iHS	SLITRK5
13	50.91	50.93	FST	ADRA1D
15	45.72	46.28	iHS	ENSBTAG00000055271, ZNF214, OR2D2, OR6A2, ENSBTAG00000051921, ENSBTAG00000027525, ENSBTAG00000049727, LOC516636, ENSBTAG00000053544, LOC100336956, ENSBTAG00000055066, U6, ENSBTAG00000051394, LOC101903126
15	46.70	47.46	iHS	ENSBTAG00000049294, CAVIN3, FAM160A2, OR52B2, OR52L1, LOC100848074, OR56A1, OR56A5, LOC613909, OR52E8, LOC523394, U6, ENSBTAG00000050950, OR52N1, ENSBTAG00000052528, LOC525863
17	70.67	71.81	FST	ENSBTAG00000050653, ENSBTAG00000053220, VPREB2, ENSBTAG00000039237, ENSBTAG00000052689, ENSBTAG00000054609, ENSBTAG00000053980, ENSBTAG00000048321, ENSBTAG00000047986, ENSBTAG00000052233, ENSBTAG00000050165, ENSBTAG00000046322, ENSBTAG00000048030, ENSBTAG00000048730, ENSBTAG00000045810, ENSBTAG00000050136, ENSBTAG00000051483, ENSBTAG00000031160, ENSBTAG00000049563, ENSBTAG00000054090, LOC100847119, ENSBTAG00000051734, ENSBTAG00000054009, LOC100847119, ENSBTAG00000048832, VPREB1, ENSBTAG00000050062 C17H22orf15, MMP11, SMARCB1, SLC2A11, MIF, GSTT4, GSTT1, GSTT2, DDT, CABIN1, bta-mir-2893, SUSD2, GUCD1
18	13.19	14.19	iHS	JPH3, BANP, ZNF469, ZFPM1, ZC3H18, CTU2, CDT1, TRAPPC2L
18	14.47	15.47	FST	SPG7, RPL13, ENSBTAG00000042786, CPNE7, ENSBTAG00000054975, DPEP1, ENSBTAG00000052192, CDK10, ZNF276, SPIRE2, TCF25, MC1R, TUBB3, DEF8, ENSBTAG00000038051, GAS8, ENSBTAG00000025283, ORC6, GPT2, ENSBTAG00000054464, ENSBTAG00000004817
18	14.93	16.09	iHS	ORC6, GPT2, U6, PHKB, ENSBTAG00000054464, ENSBTAG00000004817, ENSBTAG00000053070, ENSBTAG00000053690
18	39.08	39.96	FST	PKD1L3, ZNF821, ENSBTAG00000044454, AP1G1, ENSBTAG00000043472, PHLPP2, TAT, ZNF19, ZNF23, ENSBTAG00000053994, CMTR2, ENSBTAG00000049407, ENSBTAG00000051512
18	59.96	60.16	FST	ENSBTAG00000025023, ENSBTAG00000045985, ENSBTAG00000051883, ENSBTAG00000004925, ENSBTAG00000055311, ZNF677, ENSBTAG00000030440
20	31.35	32.35	iHS	TMEM267, CCL28, HMGCS1, SELENOP, bta-mir-12004
23	25.48	26.16	FST	ENSBTAG00000045034, BOLA-DQB, BLA-DQB, BOLA-DRB3, ENSBTAG00000048364, ENSBTAG00000038397, DSB, BTNL2, LOC504295, ENSBTAG00000026163, ENSBTAG00000050817 ENSBTAG00000048364
23	29.12	30.12	FST	OR11A12, OR11A1, LOC514434, OR12D3, OR12D21, OR12D20, LOC785431, OR5V1, LOC785162, LOC520002, ENSBTAG00000052969, LOC516273, LOC516274, LOC782301, ENSBTAG00000054760, LOC785479, ENSBTAG00000049353, LOC784681, ENSBTAG00000053092, ENSBTAG00000051060, ENSBTAG00000054399, ENSBTAG00000053389, OR10AL43, OR10AL40, ENSBTAG00000050941, LOC509280, OR2W1, LOC782475, LOC782554, ZNF311, TRIM27, U6
26	51.13	51.98	FST	ENSBTAG00000050527, JAKMIP3, ENSBTAG00000049137, DPYSL4, ENSBTAG00000050923, LRRC27, ENSBTAG00000052910, PWWP2B
27	1.12	1.81	FST	ENSBTAG00000051728, CLN8, ARHGEF10, bta-mir-10169
27	2.81	2.82	FST-iHS	U6, ENSBTAG00000049527
29	5.50	5.54	iHS	CHORDC1, TRIM77, U6, LOC616911, ENSBTAG00000052850, ENSBTAG00000052563, ENSBTAG00000053004, ENSBTAG00000048894
29	7.00	7.19	iHS	GRM5

In addition to the high number of the identified genes which were associated with olfactory receptors, 4 of the identified genes were related to gustatory receptors and the perception of tastes and flavors. Three of these genes including TAS2R39, TAS2R40, and TAS2R60 were located in the range of 106 to 106.7 Mb of chromosome BTA4. Moreover, the PKD1L3 gene was located on chromosome 18. AREG and EREG were the other important genes in selected area. These genes are involved in the activities of the growth factors and epidermal growth factors. Furthermore, the AREG gene is related to the growth and development of the breast tissue [45]. A gene, which is called CLN8, is located on chromosome BTA27 of cows. It is one of the genes that play important biological roles in the regulation of cell size, regulation of the protein catabolic processes, growth of skeletal muscles, and development of the nervous system. Moreover, it performs functional roles such as the control of social behavior and walking behavior [45]. This gene was introduced as selected gene in buffalo populations [23]. The TRIM77, PTPN22 and STAT2, RNF41, STAT6, and DDIT3 genes play a key role in the initiation of the immune response. Moreover, they participate in other biological pathways. Ghoreishifar et al. (2020) noted that the STAT6 gene played a role in the reproductive activities [46]. Furthermore, STAT6 gene has a positive role in the regulation of heat production in the body in cold conditions [45]. In addition, STAT6 and AGAP2 genes are involved in the growth and development of the mammary glands. The DDIT3 gene plays an indirect role in the production of milk by signaling the Wnt pathway. This gene was identified as a signature of selection in buffalo populations and it was suggested that the DDIT3 gene may be linked to milk production by intervening in the immune system and by triggering the immune response to mastitis [23]. The CHD8 gene is one of the genes which are involved in Wnt signaling pathway. In addition to the aforementioned functions of the DDIT3, this gene plays an important role in the development of the sense of starvation and the sensory perception of sounds. The MYO1A gene is one of the genes which are associated with the sensory perception of sounds [45].

The AGAP2 gene is one of the genes which fulfill a key role in the production of milk. It was identified as a signature of selection in the present study. AGAP2 is another gene which is involved in the growth and development of mammary alveoli. Moreover, this gene has a negative impact on the process of apoptosis and a positive effect on the JAK-STAT pathway. The JAK-STAT pathway is involved in growth, cell structure, and apoptosis [47].

AVIL, STAC3 and CYP27B1 genes were the other important genes in this analysis. The AVIL and STAC3 genes are involved in the growth and development of skeletal muscles and the development of nerve cells (HGNC Symbol). The CYP27B1 gene is involved in the immune response and bone mineral accumulation. Moreover, it plays a negative role in cell growth [45]. AGAP2 is another gene which is involved in the growth and development of alveoli of the mammary tissue. In addition, this gene has a negative effect on the process of apoptosis and a positive impact on the JAK-STAT pathway. The JAK-STAT pathway is involved in growth, cell structure, and apoptosis [47]. The AGAP2 gene prevents cell death using this pathway (HGNC Symbol).

The DNAJC14 gene is located in the 56 Mb region of bovine chromosome 5 and encodes the Heat Shock Protein (Hsp40). This gene is the co-chaperone of Hsp70 (a protein with a molecular weight of 70 kDa which is one of the most important heat shock proteins) and can have combined effects on biological functions. For instance, a gene of the same family with Hsp70 (DNAJA1 / Hsp70 combination) directly inhibits apoptosis (cell death) [48]. The CHORDC1 gene and the PTGES3 Hsp90-binding protein are the other genes which are related to heat shock, that were identified as a signature of selection in the present study.

Moreover, it is possible to mention a number of genes including ZFPM1 and ENSBTAG00000053690 genes which are involved in heart development, SNAP23 and SLITRK5 genes that are related to nervous formation and development, CHD8 and SELENOP genes which play a role in brain development, and JPH3 and GRM5 genes that are involved in learning, memory and exploratory behavior (HGNC Symbol). Furthermore, the HMGCS1, APOF and APON genes are involved in fat and cholesterol biosynthesis and the B4GALNT1 gene fulfills a role in fat storage [45]. A number of genes including ATP5F1B (involved in lung developmental metabolism), TIMELESS, SALL2 (involved in eye development), MC1R (involved in the development of coating color), and METTL3 (involved in producing responses to UV ray) are regarded to be functional and important genes. These genes were identified as the signatures of selection in the present study. Another classification is related to two genes that were melanin-related genes in the present study.

The examination of the gene networks in the DAVID and String software confirmed the existence of two noticeable gene networks which were related to olfactory traits, cell membranes and immunity. The identified QTLs from selected regions were related to milk production, milk somatic cells, fertility, ion disease, calving and growth.

Biological pathway analysis

The information on the gene groups, which were related to the identified genes using the F_ST and iHS methods, is provided in Table 2. As shown in Table 2, the number of genes in some biological pathways was high and reached 83 genes in one biological pathway. The information on a number of genes was not accessible on DAVID database. Nonetheless, 21 gene categories were determined using the information on the identified genes on DAVID. These 21 categories could be classified into three main groups. More details about these results are provided below (see Table 2). The information on the identified categories involved the number and percentage of the genes which were involved in each category, significance indices including the p value, benfroni, benjamin and FDR. As shown in the Table, most of the gene groups were clearly significant (p < 2.5E-15). This issue highlighted the existence of very strong genetic connections within each of the gene categories. The obtained information highlighted the fact that, the identified genes played roles in immune response pathways, olfactory system, G protein signaling pathway, regulation of cell proliferation, collagen catabolic process, response to lipopolysaccharide, glutathione metabolism, inflammatory response, and melanin biosynthesis. More detailed information on this section is provided in Table 2.

[Table 2 Position]

The results of the immune response pathways, the olfactory system, and the G protein signaling pathway are very important among the results of the analysis regarding the ontology of the identified genes due to their high number of genes and their high significance. They will be examined in detail below. Considering the identification of gene networks with high significance and high enrichment score (up to 9.51), it was necessary to examine the networks. Therefore, the gene networks were visually examined using the String database [49]. The extracted gene networks from the gene clusters indicated the existence of at least three very extensive networks among the genes which were identified as the signatures of selection. One of these gene networks, which included the first to the third classifications in Table 2, involved the olfaction, sensory transduction, and G-protein coupled receptor categories (Fig. 5).

Table 2

The gene categories were determined using the information on the identified genes using DAVID.
row		Term	Genes Count	gene %	P-Value	Fold Enrichment	Bonferroni	Benjamini	FDR
1	UP_KEYWORDS	Olfaction	43	18.1	2.50E-15	4.2	4.10E-13	4.10E-13	4.00E-13
2	UP_KEYWORDS	Sensory transduction	44	18.6	7.40E-15	4	1.20E-12	6.20E-13	5.90E-13
3	UP_KEYWORDS	G-protein coupled receptor	46	19.4	8.70E-13	3.4	1.50E-10	4.90E-11	4.70E-11
4	UP_KEYWORDS	Transducer	46	19.4	3.80E-12	3.2	6.40E-10	1.60E-10	1.50E-10
5	UP_KEYWORDS	Receptor	48	20.3	4.00E-09	2.5	6.70E-07	1.30E-07	1.30E-07
6	UP_KEYWORDS	Cell membrane	49	20.7	1.70E-08	2.4	2.80E-06	4.70E-07	4.50E-07
7	UP_KEYWORDS	Chemotaxis	7	3	1.50E-05	12.8	2.60E-03	3.70E-04	3.60E-04
8	UP_KEYWORDS	Cytokine	9	3.8	1.60E-04	5.8	2.60E-02	3.30E-03	3.20E-03
9	UP_KEYWORDS	MHC II	4	1.7	1.50E-03	17.1	2.20E-01	2.80E-02	2.70E-02
10	UP_KEYWORDS	Metalloprotease	5	2.1	8.40E-03	6.2	7.60E-01	1.40E-01	1.40E-01
11	UP_KEYWORDS	Transmembrane helix	80	33.8	1.60E-02	1.2	9.30E-01	2.30E-01	2.20E-01
12	UP_KEYWORDS	Transmembrane	80	33.8	1.70E-02	1.2	9.40E-01	2.30E-01	2.20E-01
13	UP_KEYWORDS	Secreted	20	8.4	1.80E-02	1.8	9.50E-01	2.30E-01	2.20E-01
14	UP_KEYWORDS	Membrane	87	36.7	3.80E-02	1.2	1.00E + 00	4.60E-01	4.40E-01
15	UP_KEYWORDS	Zinc	22	9.3	4.40E-02	1.6	1.00E + 00	4.60E-01	4.50E-01
16	UP_KEYWORDS	Inflammatory response	4	1.7	4.40E-02	5.1	1.00E + 00	4.60E-01	4.50E-01
17	UP_KEYWORDS	Signal	47	19.8	5.20E-02	1.3	1.00E + 00	5.10E-01	4.90E-01
18	UP_KEYWORDS	Protease	9	3.8	5.40E-02	2.2	1.00E + 00	5.10E-01	4.90E-01
19	UP_KEYWORDS	Tumor suppressor	3	1.3	6.80E-02	7	1.00E + 00	5.80E-01	5.50E-01
20	UP_KEYWORDS	Nuclease	4	1.7	6.90E-02	4.2	1.00E + 00	5.80E-01	5.50E-01
21	UP_KEYWORDS	Melanin biosynthesis	2	0.8	7.50E-02	25.7	1.00E + 00	6.00E-01	5.80E-01

These gene groups are generally active in a single biological pathway or a single gene network, which is called the olfactory receptor pathway. The identification of genes, which are related to the olfactory pathway, as the indicator of the genes which are associated with evolution and domestication has been reported in a number of the aforementioned studies.

[Figure 5 Position]

Figure 5. Gene network related to the olfactory pathway and sensory transmission

The other pathway, which involved most of the categories in Table 2 such as Transducer, Receptor, and Cell membrane categories, is related to membrane receptors and biological pathways for disease resistance. The figures which are associated with this network are shown in Figs. 6 and 7.

[Figure 6 Position]

Figure 6. Gene network associated with the transmembrane helix and transmembrane pathway

[Figure 7 Position]

Figure 7. Biological pathways related to immunity and disease resistance

QTLs related to the selected genomic regions

The information on the identified QTLs using the F_ST method, the QTL codes, and the cited sources is provided in the table supplementary 1. As shown in the table supplementary 1, the QTLs, which were identified using the F_ST method, were associated with milk production, milk somatic cells, fertility, ion disease, calving and growth. Considering, the provided explanation of the identified genes, there was an acceptable overlap between the identified genes and QTLs. The information on the identified QTLs using the iHS method, QTL codes, and the cited sources is provided in appendix (table supplementary 2). As shown in the table supplementary 2, the identified QTLs were related to milk production, feed efficiency, hemoglobin level, body weight, degree of obesity, growth, reproductive traits and calving. Considering the provided explanation of the identified genes, there was an acceptable overlap between the milk and fat production genes and the identified QTLs.

The identified genes can provide a better understanding of the functional roles of genes and the breeding programs. The identified genomic regions are very useful for analyzing the evolution and domestication of buffaloes. It is suggested that the information in this section be used to analyze the genetic and evolutionary structures of populations. Determining the association of genotypic with functional traits such as milk production, growth, meat quality, fertility and susceptibility to some diseases with exact and accurate phenotypic records will give more useful information

SNPs

Single Nucleotide Polymorphisms

QTL

Quantitative Trait Loci.

PCA

Principal Component Analysis.

Linkage Disequilibrium.

F_ST

Population differentiation index

HWE

Hardy-Weinberg Equilibrium

MAF

Minor Allele

Ethical Approval

Ethical review and approval was not required as no animal work was undertaken and the data were obtained from research published by Colli et al. (2018).

Competing interests

The authors declare that they have no competing interests.

Funding

This study was supported by Urmia University with fund ID: 259801. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Author Contribution

HM MGh MM JLW conceived and designed the study, MM contributed reagents/materials/analysis tools, HM MM analyzed the data and MGh MM wrote the paper.

Data Availability

The datasets supporting the results and conclusions of our study are included within the article and in the additional files. Raw data are available at Dryad information repository (doi: 10.5061 / dryad.h0cc7).

Aminafshar M, Amirinia C, Torshizi RV. Genetic diversity in buffalo population of guilan using microsatellite markers. J Anim Vet Adv. 2008;7:1499–502.
Holsinger KE, Weir BS. Genetics in geographically structured populations: defining, estimating and interpreting F ST. Nat Rev Genet. 2009;10(9):639–50.
Simianer H, Ma Y, Qanbari S. Statistical problems in livestock population genomics. in Proceedings of the 10th World Congress on Genetics Applied to Livestock Production. 2014.
Qanbari S, et al. Classic selective sweeps revealed by massive sequencing in cattle. PLoS Genet. 2014;10(2):e1004148.
Akey JM, et al. Interrogating a high-density SNP map for signatures of natural selection. Genome Res. 2002;12(12):1805–14.
Sabeti PC, et al. Detecting recent positive selection in the human genome from haplotype structure. Nature. 2002;419(6909):832–7.
Tajima F. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 1989;123(3):585–95.
Rubin C-J, et al. Whole-genome resequencing reveals loci under selection during chicken domestication. Nature. 2010;464(7288):587–91.
Nielsen R, Yang Z. Likelihood models for detecting positively selected amino acid sites and applications to the HIV-1 envelope gene. Genetics. 1998;148(3):929–36.
Grossman SR, et al. A composite of multiple signals distinguishes causal variants in regions of positive selection. Science. 2010;327(5967):883–6.
Biswas S, Akey JM. Genomic insights into positive selection. Trends Genet. 2006;22(8):437–46.
Hayes BJ, et al. A genome map of divergent artificial selection between Bos taurus dairy cattle and Bos taurus beef cattle. Anim Genet. 2009;40(2):176–84.
Yurchenko AA, et al. Scans for signatures of selection in Russian cattle breed genomes reveal new candidate genes for environmental adaptation and acclimation. Sci Rep. 2018;8(1):12984.
Igoshin A, et al. Genome-wide association study for body weight in cattle populations from Siberia. Anim Genet. 2019;50(3):250–3.
Weldenegodguad M, et al. Whole-genome sequencing of three native cattle breeds originating from the northernmost cattle farming regions. Front Genet. 2019;9:728.
Salehi R, et al. Genomic Selection Signatures in Two French and Swedish Holstein Cattle Breeds Provide Evidence for Several Potential Candidate Genes Linked to Economic Traits. Iran J Appl Anim Sci. 2023;13(4):677–84.
Biabani P et al. Detection of genetic differences between Holstein and Iranian north-west indigenous hybrid cattles using genomic data. 2022.
Brito LF, et al. Genetic diversity and signatures of selection in various goat breeds revealed by genome-wide SNP markers. BMC Genomics. 2017;18:1–20.
Shi H, et al. Whole genome sequencing revealed genetic diversity, population structure, and selective signature of Panou Tibetan sheep. BMC Genomics. 2023;24(1):50.
Dzomba E, et al. Selection signature analysis and genome-wide divergence of South African Merino breeds from their founders. Front Genet. 2023;13:932272.
Hlongwane NL, et al. Identification of Signatures of Positive Selection That Have Shaped the Genomic Landscape of South African Pig Populations. Animals. 2024;14(2):236.
Munoz M, et al. Genomic diversity, linkage disequilibrium and selection signatures in European local pig breeds assessed with a high density SNP chip. Sci Rep. 2019;9(1):13546.
Mokhber M, et al. A genome-wide scan for signatures of selection in Azeri and Khuzestani buffalo breeds. BMC Genomics. 2018;19:1–9.
Strillacci MG, et al. A genome-wide scan of copy number variants in three Iranian indigenous river buffaloes. BMC Genomics. 2021;22(1):305.
Sun T, et al. Selection signatures of Fuzhong Buffalo based on whole-genome sequences. BMC Genomics. 2020;21:1–10.
Grilz-Seger G, et al. Genome-wide homozygosity patterns and evidence for selection in a set of European and near eastern horse breeds. Genes. 2019;10(7):491.
Han H, et al. Selection signatures for local and regional adaptation in Chinese Mongolian horse breeds reveal candidate genes for hoof health. BMC Genomics. 2023;24(1):35.
Bahbahani H, et al. Positive selection footprints and haplotype distribution in the genome of dromedary camels. animal. 2024;18(3):101098.
Al Abri M et al. Assessing genetic diversity and defining signatures of positive selection on the genome of dromedary camels from the southeast of the Arabian Peninsula. Front Veterinary Sci, 2023. 10.
Almeida OAC, et al. Identification of selection signatures involved in performance traits in a paternal broiler line. BMC Genomics. 2019;20:1–20.
Asgari Z, et al. Bayes factors revealed selection signature for time to market body weight in chicken: a genome-wide association study using BayesCpi methodology. Italian J Anim Sci. 2021;20(1):1468–78.
Colli L, et al. New insights on water buffalo genomic diversity and post-domestication migration routes from medium density SNP chip data. Front Genet. 2018;9:53.
Purcell S, et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007;81(3):559–75.
Elhaik E. Principal component analyses (PCA)-based findings in population genetic studies are highly biased and must be reevaluated. Sci Rep. 2022;12(1):14683.
Aulchenko YS, et al. GenABEL: an R library for genome-wide association analysis. Bioinformatics. 2007;23(10):1294–6.
Alexander DH, Novembre J, Lange K. Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 2009;19(9):1655–64.
Gautier M, Vitalis R. rehh: an R package to detect footprints of selection in genome-wide SNP data from haplotype structure. Bioinformatics. 2012;28(8):1176–7.
Browning BL, et al. Fast two-stage phasing of large-scale sequence data. Am J Hum Genet. 2021;108(10):1880–90.
Sherman BT et al. DAVID: a web server for functional enrichment analysis and functional annotation of gene lists (2021 update). Nucleic acids research, 2022. 50(W1): pp. W216-W221.
Szklarczyk D, et al. The STRING database in 2021: customizable protein–protein networks, and functional characterization of user-uploaded gene/measurement sets. Nucleic Acids Res. 2021;49(D1):D605–12.
Mokhber M, et al. Detection of selection signatures in Azeri and Mazandrani buffalo populations by high density SNP markers. Iran J Anim Sci. 2019;50(2):89–102.
Moreno-Estrada A, et al. Signatures of selection in the human olfactory receptor OR5I1 gene. Mol Biol Evol. 2008;25(1):144–54.
Chen R, Irwin DM, Zhang Y-P. Differences in selection drive olfactory receptor genes in different directions in dogs and wolf. Mol Biol Evol. 2012;29(11):3475–84.
Groenen MA, et al. Analyses of pig genomes provide insight into porcine demography and evolution. Nature. 2012;491(7424):393–8.
Maglott D, et al. Entrez Gene: gene-centered information at NCBI. Nucleic Acids Res. 2005;33(suppl1):D54–8.
Ghoreishifar SM, et al. Genomic measures of inbreeding coefficients and genome-wide scan for runs of homozygosity islands in Iranian river buffalo, Bubalus bubalis. BMC Genet. 2020;21:1–12.
Rosenfeld RG, et al. Identification of the first patient with a confirmed mutation of the JAK-STAT system. Pediatr Nephrol. 2005;20:303–5.
Gotoh T, et al. hsp70-DnaJ chaperone pair prevents nitric oxide-and CHOP-induced apoptosis by inhibiting translocation of Bax to mitochondria. Cell Death Differ. 2004;11(4):390–402.
Szklarczyk D, et al. The STRING database in 2023: protein–protein association networks and functional enrichment analyses for any sequenced genome of interest. Nucleic Acids Res. 2023;51(D1):D638–46.

No competing interests reported.

Download PDF

Version 1

posted

You are reading this latest preprint version

Detection of Selection Signatures in some of the Water Buffaloes across the World

Status:

Version 1

Abstract

Figures

Introduction

Materials and Methods

Genomic Data

Quality Control and Data Screening

Principal Component Analysis and Population Structure

Detection of Selection Signatures

Biological Pathway Analysis

Results and Discussion

Quality control and data screening

Principal Component Analysis (PCA) and population structure

Detection of selection signatures

QTLs related to the selected genomic regions

Biological pathway analysis

QTLs related to the selected genomic regions

Conclusion

Abbreviations

Declarations

Ethical Approval

Competing interests

Funding

Author Contribution

Data Availability

References

Additional Declarations

Supplementary Files

Status:

Version 1