Comparative Genomic Characterization of Insulin-Like Growth Factor Binding Proteins in Cattle and Buffalo.

The somatotropic axis consists of genes associated with economic traits like muscle growth and carcass traits in livestock. Insulin-like growth factor binding proteins (IGFBPs) are the major proteins that play a vital role in the somatotropic axis. The present study performed a genome-wide characterization of IGFBP genes in cattle. Genomic sequences of the IGFBP gene family for different mammals (cattle, buffalo, goat, and sheep) were recovered from the NCBI database. Sequence analyses were performed to investigate cattle's genomic variations in the IGFBP gene family. Phylogenetic analysis, gene structure, motif patterns, and conserved domain analysis (CDA) of the IGFBP family revealed the evolutionarily conserved nature of the IGFBP genes in cattle and other studied species. Physicochemical properties of IGFBP proteins in cattle revealed that most of these proteins are unstable, hydrophilic, thermostable, and acidic. Comparative amino acid analysis revealed variations in all protein sequences with single indels in IGFBP3 and IGFBP6. Mutation analysis revealed only one nonsynonymous mutation D212 > E in the IGFBP6 protein of cattle. A total of 245 nuclear hormone receptor (NHRs) sites were detected, including 139 direct repeats (DR), 65 everted repeats (ER), and 41 inverted repeats (IR). Out of 133 transcription factors (TFs), 10 TFs (AHR, AHRARNT, AP4, CMYB, E47, EGR2, GATA, SP1, and SRF) had differential distribution (P value < 0.05) of putative transcriptional binding sites (TFBS) in cattle compared to buffalo. Synteny analysis revealed the conserved nature of genes between cattle and buffalo. Two gene pairs (IGFBP1/IGFBP3 and IGFBP2/IGFBP5) showed tandem duplication events in cattle and buffalo. This study highlights the functional importance of genomic variation in IGFBP genes and necessitates further investigations better to understand the role and mechanisms of IGFBPs in bovines.

1. Introduction

IGFBPs are cysteine-rich proteins with high essential amino acid sequence similarity. Almost six mammalian types of IGFBPs commonly exist, ranging from IGFPB1 to 6. They can bind insulin-like growth factor (IGF) growth with a very high affinity [1]. IGF-independent activities of IGFBPs significantly affect biological procedures, including angiogenesis or cell proliferation [1, 2]. The expression of IGFBP is controlled by various physiological conditions such as exercise, nutrition, and pregnancy. IGFBPs are used as a biomarker for improving husbandry conditions, health status, phenotype selection, or disease analysis in farm animals [3].

IGFs are the most important polypeptides in mammals that control body growth, metabolism, mitosis, and cell differentiation [4, 5]. IGFs have a galactopoietic role and stimulate the development and growth of the mammary glands [6]. IGFBPs are a member of a family of six homologous proteins that always bind with IGFs. They control most of the biological activities (metabolism and growth). Based on amino acid similarity analysis, IGFBP1 is more closely related to IGFBP2 than to IGFBP3, which in turn is more closely associated with IGFBP5 [1]. IGFBP3 induces cell proliferation or apoptosis (programmed cell death) in various cells [7]. Lower levels of IGFBP2 in estrogen-dominant bovine follicles are not attributable to enhanced proteolysis; lower levels of IGFBP4 and 5 are most likely to exhibit protease activity [8]. IGFBP4 protease plays a significant role in developing dominant ovarian follicles [9]. IGFBP6 levels, on the other hand, were shown to rise with age [10].

The IGFBP3 is a direct growth inhibitor and a regulator of IGF bioactivity in extravascular tissues [11]. IGFBP3 was the primary component of six IGFBPs in mammalian circulation [12]. In plasma, it is associated with more than 90% of IGFs [13]. The IGFBP3 gene is present on chromosome 4 in cattle [14] and is 8.9 kb in length, having five exons [15]. The IGFBP3 is involved in various body functions, including metabolism, immunity, growth, energy balance, and reproduction. Due to the primary role of IGFBP3 in animal growth and maturity, the IGFBP3 gene is considered a candidate gene and a good marker for production and growth traits [16].

The development of the IGFBP gene family has been attributed to repetitive chromosomal duplications. According to thorough sequence analysis, sequence-based phylogenies, and chromosomal information, the ancestral chordate IGFBP gene experienced a localized gene duplication, forming a gene pair close to a HOX cluster [17]. The IGFBP gene family follows the six IGFBP types in today's placental mammals due to 2 basal vertebrate tetraploidization. Many IGFBP genes have survived despite that the gene family's ancient expansion strongly implies that each gene gained different and significant functions early in mammal evolution [17].

The genomic characterization of the IGFBP gene family is necessary to provide insights into genetic variation in different IGFBP genes and their comparative physiological role in production performance in cattle. Moreover, identifying transcription binding sites and nuclear receptors can help better comprehend the regulation of IGFBP genes in cattle. Therefore, the present study was planned for the genomic characterization of the IGFBPs in cattle to elucidate nucleotide sequence variations.

2. Materials and Methods

2.1. Identification of IGFBP Genes in Cattle

The cattle (ARS-UCD1.2), buffalo (UOA_WB_1), sheep (Oar Rambouillet_v1.0), and goat (ARS1) genomes, proteomes, and annotation data were retrieved from the National Center for Biotechnology Information (NCBI) genome database [18] aiming at identifying IGFBP genes in cattle. The GenBank accession number of each member of the IGFBP gene family of every species is listed in supplementary Table S1. To validate all potential IGFBP genes in cattle at a genome-wide level, both the basic local alignment search tool (BLAST) and the hidden Markov model (HMM) searches were carried out. The buffalo (Bubalus bubalis), goat (Capra hircus), and sheep (Ovis aries) IGFBP (insulin-like growth factor binding protein) sequences were also validated from the UniProt database [19].

The sequences were submitted to BLASTp, with a threshold of e value = 10−5, using the BLOSUM62 matrix, which had a word size of six, a gap cost of eleven, an extension of one, and an adjustment for the conditional compositional score matrix. The cattle dataset was also searched with the HMMER [20] program using the HMM profile of the IGFBP domain (PF00219) from the Pfam database with an e value threshold of 1.0 × e−5. Duplicated sequences were excised after retrieving the necessary protein sequences to avoid ambiguity. These nonredundant sequences were examined using the Simple Modular Architecture Research Tool (SMART) to validate the presence of domains in the IGFBP protein sequences [21]. The conserved protein domains in the cattle IGFBP family were compared with the NCBI's conserved domain database (CDD).

2.2. Phylogenetic Analysis

ClustalW (protein sequence alignment tool) was used to align the amino acid sequences of IGFBP from cattle, buffalo, goats, and sheep. The aligned sequences were utilized to create a neighbour-joining (NJ) phylogenetic tree of the IGFBP gene family by the MEGA7 program, using a Poisson model with pairwise deletion and a bootstrap value of 1000 replicates.

2.3. Structural Feature Analysis

The MEME suite [22] was used to analyze the conserved motifs in the dataset. The IGFBP protein sequences were provided in FASTA format as a query, and a site distribution was chosen for each sequence with one occurrence of each site. The minimum and maximum widths of motifs were 6 and 50, respectively. The number of motifs was set as 10. All CDs and genomic sequences were imported into the Gene Structure Display Server 2.0. (GSDS) [23]. Then, the final gene structure was displayed and visualized using the cattle genome annotation file in general feature format (GFF) using in-house scripts in the TBtools software.

2.4. Characterization of Physicochemical Properties of IGFBP Proteins

The ProtParam was used to describe physical and chemical characteristics of IGFBP proteins, such as molecular weight (MW), amino acid number, isoelectric point (pI), aliphatic index (AI), instability index (II), and grand average of hydropathicity (GRAVY) as described, previously [24, 25].

2.5. Multiple Sequence Alignment

In order to identify indels and illustrate sequence variations in the IGFBP protein sequence, the Multiple Align Show [26] tool was used to perform multiple sequence alignment of the IGFBP protein sequences.

2.6. Mutation Analysis

The mutations observed in the protein sequences of four species were further analyzed through different online tools (PolyPhen-2, MUpro, PROVEAN, I-Mutant, PhD-SNP, SIFT, SNAP2, PredictSNP, Meta-SNP, and SNAP) to look for their impacts on protein structure and functions. Further, the 3-dimension structures of IGFBP6 genes of cattle, buffalo, goat, and sheep were predicted through MODELLER and the quality of the structures was checked in SAVES server through ERRAT, Verify 3D, PROVE, PROCHECK, and Whatceck. The Ramachandran plot showed mostly residues in allowed regions. No residue was present in the disallowed region. To calculate the RMSD value, the mutated cattle IGFBP6 protein structure was aligned against protein structures of IGFBP6 genes of buffalo, goat, and sheep in Moe 2021 (https://www.chemcomp.com/Products.htm) software to know the structural similarities and deviations.

2.7. Detection of TF Binding Sites

Detection of possible signals for a putative transcription binding factor was accomplished using the Promoter 2.0 prediction server [27]. As input, cattle and buffalo genomic sequences were provided in FASTA format. The site with a higher score (>1.0) is a probable predicted site. TFBIND [28] was used to find the sequence of 100 bp upstream from the increased likelihood predicted location, which was obtained and entered into the program. The weight matrix of the transcription factor database TRANSFAC R.3.4 was used to identify transcription factor binding sites.

2.8. Nuclear Hormone Receptor Site Identification

Further, for the prediction of nuclear hormone receptor binding sites, the NHR scan [29] was employed using genomic sequences in FASTA format. The cumulative probability of entering match states was 0.01 using the NHR scan.

2.9. Synteny Analysis

The whole genomes of cattle and buffalo were blasted to each other to identify collinear genes. TBtools was used to construct the dual-synteny plot, which was created by providing annotation files for both genomes, including information about collinear genes and chromosomal IDs.

Chromosomal locations of buffalo IGFBP genes were acquired from their genome resources. A genome annotation file with a general feature format (GFF) was given as an input to the MCScanX program [30] to map the physical locations of genes on chromosomes and then visualized in TB tools. Furthermore, the buffalo and cattle dual-synteny plots were plotted for IGFBP gene collinearity. Additionally, pairwise alignment of homologous gene pairs of IGFBP genes using MEGA7 [31] with the MUSCLE algorithm was used to assess the occurrences of duplications for the buffalo IGFBP gene family. The DnaSP v6.0 [32] software was also used to calculate pairwise synonymous substitutions per synonymous site (Ks) and nonsynonymous substitutions per nonsynonymous site (Ka) adjusted for multiple hits.

3. Results

3.1. Phylogenetic Analysis

The molecular phylogenetic study of representative species indicated that all of the IGFBP gene sequences were classified into seven types: IGFBP1, IGFBP2, IGFBP3, IGFBP4, IGFBP5, IGFBP6, and IGFBP7 (Figure 1). All seven genes were grouped into three sister clades. From top to bottom, clade 1 included IGFBP2, -4, and -1; clade 2 included IGFBP6, -3, and -5; and clade 3 included IGFBP7. Overall, phylogenetic analysis of the IGFBP gene family revealed that the Bos taurus is more closely related to Bubalus bubalis in all IGFBP gene family groups except IGFBP7 where it was distantly related. Similarly, Capra hircus is more closely associated with Ovis aries except for IGFBP2 and IGFBP7 where it was distantly related (Figure 1).

Figure 1

Molecular phylogenetic analysis of the IGFBP gene family (IGFBP1 to IGFBP7) in Bos taurus, Bubalus bubalis, Capra hircus, and Ovis aries.

3.2. Structural Categorization of the IGFBP Gene Family

The investigation of gene organization, motif patterns, and conserved domains in the IGFBP gene family in four studied species was also carried out to undertake structural characterization of the IGFBP gene family while considering their evolutionary connections (Figures 2(a)–2(d)). The top ten MEME-conserved motifs in the IGFBP genes were analyzed (Table 1). The MEME1 motif with a 41 amino acid sequence length was annotated as the Thyroglobulin_1 domain, while the MEME2 motif with a 50 amino acid sequence length was annotated as the IB (insulin growth factor-binding protein homologs) domain after Pfam search. Further, conserved domain analysis for all genes revealed Ig superfamily and Kazal domains (Figure 2(c)). Identified conserved domains were further confirmed using the CDD BLAST in NCBI. A gene structural analysis revealed that cattle IGFBP genes in the same group had a similar number of introns and exons, even though the untranslated regions (UTRs) upstream and downstream of the gene structure differed significantly between them (Figure 2(d)). However, different IGFBP genes displayed a differential distribution of introns and exons in their coding regions.

Figure 2

Phylogenetic relationships, gene structure, motif patterns, and conserved domain regions of the IGFBP gene family in Bos taurus, Bubalus bubalis, Capra hircus, and Ovis aries. (a) Phylogenetic relationship of 7 IGFBPs. (b) Motif pattern. (c) Conserved domain regions of the IGFBP gene family. (d) Gene structure of IGFBPs. Ten putative motifs are indicated in different colored boxes.

Table 1

Top ten differentially conserved motifs detected in the IGFBP (IGFBP1, IGFBP2, IGFBP3, IGFBP4, IGFBP5, IGFBP6, and IGFBP7) gene family.

Motif	Protein sequence	Length	Pfam domain
MEME1	YALYJPNCDKKGFYKKKQCKPSLGGKRGLCWCVDKYGGKLP	41	Thyroglobulin_1
MEME2	PGCGCCATCALREGQPCGVYTPRCAQGLRCYPPPGEEKPLHALLHGRGVC	50	IB
MEME3	SESRQETEQGPCRRELEKVLZRLKAEQLR	29	—
MEME4	TAGAGEVVRCEPCDEEALARCPPPPGSPP	29	—
MEME5	SEKIDGDPZCHTFDN	15	—
MEME6	PPLLLALLLAA	11	—
MEME7	KVFREKVTPIHVSMGSGGKKH	21	—
MEME8	ESTESGEIEEN	11	—
MEME9	RASADYVLLAEQLAA	15	—
MEME10	VERGSDEH	8	—

3.3. Physicochemical Properties of the IGFBPs

The physicochemical characteristics of the IGFBP gene in cattle were evaluated in terms of their location on the chromosome, exon count, molecular weight (Da), number of amino acids (A.A) in each peptide, aliphatic index (AI), isoelectric point (pI), instability index (II), and grand average of hydropathicity index (GRAVY) (Table 2). IGFBP1 and IGFBP3 were discovered on chromosome 4, IGFBP2 and IGFBP5 on chromosome 2, IGFBP4 on chromosome 19, IGFBP6 on chromosome 5, and IGFBP7 on chromosome 6 in the region between 250 kb, which contains a variable number of exons and a variable number of length of the gene with amino acid residues. The IGFBPs had a molecular weight ranging from 24 to 34 kDa. The IGFBPs in cattle were unstable but thermostable proteins, as shown by the presence of values greater than 60 for the aliphatic index of all IGFBPs. The pI values also indicated that the IGFBP1 protein is a slightly acidic peptide but the other IGFBP proteins (IGFBP2, IGFBP3, IGFBP4, IGFBP5, IGFBP6, and IGFBP7) are basic peptides in nature. The IGFBPs in cattle were unstable except for IGFBP2, which showed an instability index value smaller than 40. All proteins were thermostable due to higher values of the aliphatic index. Lower GRAVY values suggested cattle's hydrophilic nature of IGFBP proteins (Table 2).

Table 2

Physicochemical characteristics of the IGFBP gene family in cattle.

Gene	Chromosome	Exon count	MW (Da)	A.A	pI	AI	II	GRAVY
IGFBP1	4	4	28794.05	263	6.33	75.02	54.07	−0.295
IGFBP2	2	4	34014.92	317	7.13	72.02	36.20	−0.521
IGFBP3	4	5	31570.09	291	9.03	56.05	47.67	−0.636
IGFBP4	19	4	27890.00	258	7.10	70.00	46.51	−0.483
IGFBP5	2	4	30313.96	271	8.72	61.22	51.81	−0.632
IGFBP6	5	5	24966.92	237	8.73	51.10	63.20	−0.781
IGFBP7	6	5	29078.28	282	8.25	77.84	46.77	−0.193

MW: molecular weight in Daltons; A.A: number of amino acids; pI: isoelectric point; AI: aliphatic index; II: instability index; GRAVY: grand average of hydropathicity index.

3.4. Comparative Amino Acid Analysis of IGFBP

Comparative amino acid analysis of targeted species revealed 2 short indels in the IGFBP gene family, including a single indel in IGFBP3 and a single indel in IGFBP6 (Figure S1A–S1G).

In the IGFBP1, a single amino acid variation was observed at positions L40 > M and I249 > V in Bos taurus, Q19 > H in Bubalus bubalis, and L14 > P and E226 > D in Ovis aries. Comparatively, amino acids A, R, and K were observed at positions 146, 166, and 254, respectively, in Bos taurus and Bubalus bubalis, while amino acids T, K, and Q were observed at the same positions, respectively, in Capra hircus and Ovis aries (Figure S1A).

In the IGFBP2, a single amino acid variation was observed at positions A58 > S and H145 > Q in Bubalus bubalis, E130 > V, M233 > V, and E304 > G in Capra hircus, and L150 > Q, S175 > F, and G176 > R in Ovis aries (Figure S1B).

Only one indel was observed in IGFBP3 in Bos taurus, at position 132 > 133. A single amino acid variation was observed at position A17 > T in Bos taurus, A34 > T in Bubalus bubalis, A61 > T in Capra hircus, and E79 > V and L232 > P in Ovis aries. Comparatively, amino acids Q, R, S, S, S, H, V, F, P, R, H, M, G, G, D, Y, and M were observed at positions 107, 126, 134, 136, 158, 160, 164, 167, 169, 208, 225, 231, 236, 276, 278, 288, and 290, respectively, in Bos taurus and Bubalus bubalis, while amino acids H, S, P, P, T, R, D, S, L H, D, T, A, F, S, L, and T were observed at the same positions in Capra hircus and Ovis aries. Additionally, amino acid M was observed at position 32 in Bos taurus and Bubalus bubalis, while amino acids A and V were observed at the same position in Capra hircus and Ovis aries, respectively. Similarly, amino acids L and A were observed at positions 156 and 170, respectively, in Capra hircus and Ovis aries. In contrast, amino acids G and V were observed at position 156 in Bos taurus and Capra hircus, respectively, and amino acids I and V were regarded at position 170 in Bos taurus and Bubalus bubalis, respectively (Figure S1C).

In IGFBP4, single amino acid variation was observed at positions V7 > M in Bos taurus, A161 > C and P162 > T in Bubalus bubalis, and A8 > T in Ovis aries. Two amino acid variations were seen at position 166 in the targeted species. Comparatively, amino acid A was observed at position 166 in Bos taurus and Bubalus bubalis, while amino acid V was regarded at the same position in Capra hircus and Ovis aries (Figure S1D).

In IGFBP5, single amino acid variation was observed at positions P16 > S and E185 > K in Bos taurus and G258 > R in Capra hircus. Comparatively, amino acid M was observed at position 182 in Bos taurus and Bubalus bubalis, while amino acid L was observed at the same position in Capra hircus and Ovis aries (Figure S1E).

In IGFBP6, one indel was observed at position 46 in Capra hircus and Ovis aries. Single amino acid variation was observed at positions D212 > E in Bos taurus, N147 > S in Bubalus bubalis, G103 > V in Capra hircus, and Q97 > R, G190 > D, and S235 > G in Ovis aries. Comparatively, amino acid S was observed at positions 144 and 146 in Bos taurus and Bubalus bubalis, while amino acid V and P were observed at the same positions in Capra hircus and Ovis aries, respectively (Figure S1F).

In IGFBP7, single amino acid variation was observed at positions P5 > L, L21 > P, R78 > G, V229 > A, and Q256 > R in Ovis aries. Comparatively, amino acid R was observed at position 95 in Bos taurus and Ovis aries, while amino acid K was observed at the same position in Bubalus bubalis and Capra hircus (Figure S1G).

3.5. Mutation Analysis

After the comparative amino acid analysis, the mutations observed were further analyzed using online tools to predict their effects on protein structures and functions. The results differed for different tools (Table S2). For all the mutations, the overall impact was neutral (synonymous mutations) except for position D212E in IGFBP6 protein in cattle, which was found damaging (nonsynonymous).

Further, the mutated cattle IGFBP6 protein structure was superimposed against IGFBP6 protein structures of buffalo, goat, and sheep to know the structural similarities and differences (Figure S2A-S2C) and the root means square deviation (RMSD) values were calculated (Table 3). The RMSD values were less than 2 Å for all superimposed structures. “Cattle and goat” had zero RMSD values showing the identical structure, following the “cattle and buffalo” (0.073) and “cattle and sheep” (0.171) (Table 3).

Table 3

Effects of mutations on the IGFBP6 protein structure of cattle aligned against buffalo, goat, and sheep IGFBP6 protein structures.

Structure A	Structure B	RMSD value
IGFBP6 cow (blue)	IGFBP6 buffalo (yellow)	0.073
IGFBP6 cow (blue)	IGFBP6 goat (cyan)	0.000
IGFBP6 cow (blue)	IGFBP6 sheep (green)	0.171

3.6. Transcription Factor Binding Sites

The putative transcription factor binding sites (TFBS) were observed for IGFBP genes in cattle and buffalo. 133 transcriptional factors (TFs) were screened for differential distribution of TFBS, out of which only 10 TFs (AHR, AHRARNT, AP4, CMYB, E47, EGR2, GATA, SP1, and SRF) were found to have the high differential distribution of TFBS. Both cattle and buffalo had a variable number of binding sites for most TFs (Table S3).

3.7. NHR Sites in IGFBP

The pattern of NHR (nuclear hormone receptor) sites in the IGFBP gene family in cattle was investigated using the genomic sequence of the Bos taurus. 245 NHR sites were identified in the cattle IGFBP gene family (Figure S3). Moreover, the identified NHRs in IGFBP1, IGFBP2, IGFBP3, IGFBP4, IGFBP5, IGFBP6, and IGFBP7 were 6, 50, 12, 26, 28, 6, and 117, respectively. In total, 139 direct repeats (DR) and 65 everted repeats (ER) were found in the cattle IGFBP genes which are prominently used by type II receptors (RXR) and some type III receptors (orphan receptors) can also use DR. The DR distributed in IGFBP1, IGFBP2, IGFBP3, IGFBP4, IGFBP5, IGFBP6, and IGFBP7 were 3, 31, 8, 15, 17, 4, respectively, and 61 and ER were 3, 13, 2, 3, 6, 0, and 38, respectively. A total of 41 inverted repeats (IR) were observed in different IGFBP genes primarily used as the hormonal response element (HRE) important for steroid receptors. The IR distributed in IGFBP1, IGFBP2, IGFBP3, IGFBP4, IGFBP5, IGFBP6, and IGFBP7 were 0, 6, 2, 8, 5, 2, and 18, respectively.

3.8. Synteny Analysis

Collinearity analysis showed that IGFBP genes were randomly distributed over 5 chromosomes in cattle and buffalo (Figure 3). In cattle, IGFBP genes were present on chromosomes 2, 4, 5, 6, and 19, whereas these genes were located on chromosomes 2, 3, 4, 7, and 8 in buffalo. Most genes were present close to the centromere of chromosomes in cattle.

Figure 3

Synteny plot between cattle and buffalo genomes.

Gene duplication events were observed in IGFBP family members of cattle and buffalo using Tb tools (Figures 4(a) and 4(b)). No segmental duplication was observed in both species. Tandem duplication was observed between homologous gene pairs IGFBP1/IGFBP3 and IGFBP2/IGFBP5 for both species at chromosome positions 4 and 2 for cattle and 8 and 2 for buffalo, respectively. Further, the ratios of nonsynonymous substitutions per nonsynonymous site (Ka) to synonymous substitutions per synonymous site (Ks) were calculated for these duplication events. The results indicated that duplicated gene pairs IGFBP1/IGFBP3 and IGFBP2/IGFBP5 had 0.81 and 0.86 Ka/Ks ratios in cattle and 0.79 and 0.73 in buffalo, respectively (Table 4).

Figure 4

Gene duplication events in the IGFBP gene family of (a) cattle and (b) buffalo.

Table 4

Analysis of the Ka/Ks ratios for each pair of the duplicated IGFBP gene pair in cattle and buffalo.

Cow
Gene pairs	Chr	Duplication	Ka	Ks	Ka/Ks
IGFBP1/IGFBP3	4/4	TD	0.7143	0.8809	0.81
IGFBP2/IGFBP5	2/2	TD	0.6230	0.7288	0.86
Buffalo
Gene pairs	Chr	Duplication	Ka	Ks	Ka/Ks
IGFBP1/IGFBP3	8/8	TD	0.6848	0.8641	0.79
IGFBP2/IGFBP5	2/2	TD	0.6021	0.8202	0.73

4. Discussion

Recent developments in genomic sequencing technology, particularly next-generation sequencing, resulted in the availability of sequenced genomes for many animal species, opening new avenues for understanding the genomic architecture at the molecular level [33]. Comparative genomics allows for identifying novel genes and the functional components underlying the unique traits of different species [34].

4.1. Evolutionary and Phylogenetic Analysis

Understanding the genetics and evolutionary processes of physiologically crucial genes, such as the IGFBP gene family in mammals, is necessary for understanding the regulatory mechanisms of these genes. Previously, the IGFBP gene family has been reported to contain six members (IGFBP1–IGFBP6) in vertebrates [1, 17]. Our molecular phylogenetic analysis of the IGFBP gene family revealed that all the four representative species contained seven genes (IGFBP1–IGFBP7) (Figure 1). The Bos taurus was closely related to Bubalus bubalis, whereas Capra hircus and Ovis aries shared higher sequence similarities. Our results are also supported by the study performed in Capra hircus, where Ovis aries and Capra hircus were closely related evolutionary, whereas Bos taurus was distantly related [35].

4.2. Structure of the IGFBP Gene Family

The top ten MEME-conserved motifs were observed in the IGFBP genes in the present study. The six human IGFBPs ranged in size between 240 and 328 residues and shared a common structural organization with two conserved domains separated by a variable central region. The N-terminal domain contains the primary IGF-binding (IB) site, and the C-terminal domain is a thyroglobulin type-1 domain [11, 17]. Our study annotated the MEME1 motif with a 41 amino acid sequence length as the thyroglobulin_1 domain. Previously, the thyroglobulin_1 domain is present in several protein families, including IGFBP6, a member of IGFBPs [36]. Thyroglobulin-1 repeats exhibit inhibitory activity against cysteine proteases and show greater selectivity in their interactions with target proteases. The MEME2 motif with a 50 amino acid sequence length was annotated as IB (insulin growth factor-binding protein homologs) domain after Pfam search (Table 1). The presence of thyroglobulin_1 and IB domains in C-terminal and N-terminal regions, respectively, suggested that the IGFBPs across all the studied species remained conserved throughout the evolution. Further, similar intron and exon structures of IGFBP genes across all studied species supported the conserved structural organization of IGFBPs. The upstream and downstream UTR structure considerably varied in different IGFBP genes mainly attributed to the absence or presence of retroposon elements.

4.3. IGFBP Physicochemical Properties

The physicochemical properties of the IGFBP gene family in cattle revealed that their molecular weight ranged from 24 to 34 kDa. The pI values indicated that all proteins are basic except IGFBP1, slightly acidic. The instability index (II) reflects the stability of proteins within the test tube. The II values lower than 40 indicate stable proteins [37]. The IGFBP peptides in cattle were unstable except for IGFBP2, which showed an instability index value smaller than 40. Higher values of AI reveal the thermostability of globular proteins [38]. The AI values of IGFBP1, IGFBP2, and IGFBP4 were more than 65, indicating their higher thermostability than IGFBP3, IGFBP5, and IGFBP6. All proteins were found thermostable due to higher values of the aliphatic index. GRAVY values tell about the hydropathic character of proteins [39]. In our study, negative GRAVY values suggested that cattle's IGFBP proteins are hydrophilic.

4.4. Comparative Amino Acid Sequence Analysis

Comparative genomics is a large-scale, integrative method that analyzes two or more genomes to determine their similarities and differences and investigate the biology of the individual genomes. To acquire diverse views on the organisms, comparative studies may be carried out at different levels of the genomes [40].

Comparative amino acid analysis of cattle with buffalo, goat, and sheep revealed amino acid variations in all IGFBP protein sequences (Figure S1). Overall, all six proteins were well conserved in all four analyzed species with a size range of 263 to 317 amino acid residues. The IGFBP superfamily was also found conserved in humans with a size range of 240 to 328 amino acid residues [17]. Deletions were observed in IGFBP3 in cattle at position 132 > 133 and in IGFBP6 in goats and sheep at position 46. Insertions, deletions, and mutations have played an important role in IGFBP family members' divergence from their progenitor [41]. Most of the variations were observed in the IGFBP3 gene. IGFBP3 is the most common IGF carrier protein in bovine serum, binding more than 95 percent of the IGF-I and II in circulation [42]. The nucleotide sequencing study revealed similarity percentages of 88.54, 89.63, and 95.0 in IGFBP3 gene sequences between “sheep and cattle,” “sheep and buffalo,” and “cattle and buffalo,”, respectively [43]. In another study, sheep's IGFBP3 gene nucleotide and amino acid sequence were compared with cattle and buffalo and 90 and 93 percent similarities were observed, respectively [44].

4.5. Mutation Analysis

Mutations can improve the overall fitness or can damage the structural conformation resulting in altered functions or being neutral without any change [45]. Synonymous mutations are those where no change in amino acid sequence occurs due to nucleotide change, while nonsynonymous mutations change the protein sequence by altering the amino acid. In our study, single amino acid variations observed by aligning the sequences of all target species were further analyzed to check for their effects on the protein structure and functions. Only one nonsynonymous mutation was observed at position D212E in IGFBP6 protein in cattle using different online tools (Table S1). IGFBP6 protein sequence alignment comparison of human with rat, mouse bovine, pig, and zebrafish revealed 70–85% sequence identity [46]. The same amino acid variation was also recognized in that comparison only for cattle. Mutations in the protein structure can cause problems in protein catalytic activity, stability, and interaction with other molecules [45]. The IGFBP6 role has been determined in folliculogenesis [47], carcass traits [48], and 305 days in milk yield [49] in cattle. So, a nonsynonymous mutation in IGFBP6 protein may affect these functions. The altered expression of the IGFBP6 gene in different tissues of abnormal cloned cattle was linked with increased birth weight and other abnormalities [50].

RMSD values quantify the structural similarities between two or more proteins. The present study showed lower RMSD values < 2 Å, which indicated high structural similarity in the IGFBP6 protein of cattle with buffalo and sheep, and a zero RMSD value for “cattle and goat” indicated that the structures are identical in conformation [51]. RMSD is also one of the methods to quantify the sequence alignment and evolutionary similarities between the two proteins [52, 53].

4.6. Putative Transcription Factor Binding Sites

A total of 133 TFs were screened for differential distribution of putative TFBS, out of which 10 TFs were found to have a differential distribution of binding sites in cattle and buffalo (Table S3). Aryl hydrocarbon receptor (AhRs) binding factor sites were found variable in cattle and buffalo. A study indicated that AhR can suppress the IGF pathway via upregulation of IGFBP1, resulting in decreased IGFs' bioavailability [54]. Activator protein 4 (AP4) is also an important TF and plays an important role in cell proliferation, cell growth, apoptosis, and metabolism [55]. AP4 binds to the E-box motifs at the IGFBP-2 promoter in luciferase reporter assays and serves as a transactivator [56]. C-myb is a member of the myb protein family involved in cell growth, differentiation, and apoptosis and acts as a transcriptional transactivator [57]. C-myb is highly expressed in immature hematopoietic cells and is known as an oncogene. C-myb was found to regulate the function of the IGFBP3 gene in leukaemia cells and the IGFBP5 gene in neuroblastoma cells [58, 59]. Transcription factor E47 is involved in gene regulation of muscle tissue differentiation by interacting with MyoD [60]. Specificity protein 1 (Sp1) is an important binding factor and is observed to control cell growth and its overexpression leads to tumour formation [61]. In rats, the TATA-less promoter for the IGFBP2 gene requires three clustered Sp1 sites for effective transcription [62]. The serum response factor (SRF) is related to the MADS (MCM1, Agamous, Deficiens, and SRF) family of proteins and plays an important role in skeletal muscle growth in adult mammals [63].

4.7. NHR Patterns in IGFBP

Understanding the complex biochemical systems that regulate gene transcription is essential for understanding the information flow from gene to protein and, as a result, the cell's dynamics [64]. Nuclear receptors regulate transcription by attaching DNA sequences in target genes to hormone response elements (HREs). These elements are found in regulatory sequences usually found in the target gene's 5-flanking region. Although HREs are frequently located near the main promoter, many kilobases upstream of the transcriptional start point can also be found in enhancer regions [65]. A single NHR frequently affects many genes, and various NHRs may compete for the same target sites, resulting in target gene networks that overlap [66]. NHR can also suppress the gene expression due to competition for the same target site or binding with negative HREs [65]. The pattern of NHR sites in the IGFBP gene family in cattle was investigated using the genomic sequence of the Bos taurus. 245 NHR sites were found in the cattle IGFBP gene family (Figure S2). In a total of 139 direct repeats (DR), 65 everted repeats (ER) and 41 inverted repeats (IR) were identified in the cattle IGFBP genes.

4.8. Synteny Analysis

Synteny blocks are chromosomal segments shared by two genomes with the same order of homologous genes originating from a common ancestor site [67]. The capability to examine evolutionary mechanisms that lead to diversification of a chromosomal number and shape in numerous lineages throughout the tree of life has been made possible by comparing genomic synteny between and within species [68, 69]. Collinearity analysis showed that IGFBP genes were randomly distributed over 5 chromosomes in cattle and buffalo showing the collinear relationship (Figure 3). Comparative genomic studies between cattle and buffalo revealed up to 97% similarities [70, 71]. Expansion of the genome happened during evolution through gene duplication events, leading to the increased genome size of organisms as indicated by the “2R hypothesis” [72, 73]. All the duplicated genomes did not get fixed, only 5 to 10% got fixed, and others were lost in the evolutionary process [74, 75]. In cattle and buffalo, predominantly tandem duplication events were observed, revealing the role of tandem duplication in the expansion of the IGFBP gene family. Our results are also supported by Liu et al. [76], who suggested that 75–90% of segmental duplications are organized into local tandem duplication clusters in the cattle genome. These genomic duplication events that happened during evolution helped in adaption and speciation. The duplicated genes under positive selection pressure show >1 Ka/Ks ratios, whereas duplicated genes under purifying pressure show <1 Ka/Ks ratios [77]. In our results, duplicated gene pairs in cattle and buffalo showed <1 Ka/Ks ratios, indicating purifying pressure for these genomic duplications. Due to two basic vertebrate tetraploidization, the IGFBP gene family follows the six IGFBP types found in today's placental animals. The fact that numerous IGFBP genes have survived despite the gene family's ancient expansion clearly suggested that each gene acquired distinct and important tasks early in mammalian evolution [17].

5. Conclusion

The present study concluded that the IGFBP gene family remained well conserved throughout its evolution in cattle and buffalo. Mutation analysis revealed one nonsynonymous mutation at position D212 > E in IGFBP6 in cattle which may affect important functions like folliculogenesis, growth and development, and lactation. Differential distribution of NHRs and TFBS in cattle and buffalo indicated the variation of putative regulation of these genes. Gene duplications revealed the role of tandem duplication in the expansion of the IGFBP gene family in bovines, and these genome duplication events helped in adaption and speciation.