Association Analysis of Polymorphism of FGFR1 and EBP41L5 Genes with Kidding Performance of Goats

In the development process of human civilization, goats provided humans with the main means of living such as meat, milk, skin and hair, and played an important role in human agricultural civilization and economic development. The kidding traits of goats were closely related to the reproduction speed and scale of their populations. They were an important factor that affects the economic benefits of goat breeding and were also important economic traits of goat breeding (Zheng et al., 2021). However, due to the low heritability of kidding traits, it was relatively slow to use traditional breeding methods to get more obvious genetic improvement of kidding traits, while molecular breeding technology can quickly increase the kidding rate of goats by controlling the major genes of kidding traits. Many scientists have conducted a lot of research on the relationship between single nucleotide polymorphisms and goat kidding traits, hoping to find genes that can be used to improve goat kidding traits (Zhu et al., 2021). Fibroblast growth factor receptor 1 (FGFR1) belongs to the fibroblast growth factor receptors (FGFRs) family, which includes a total of 4 receptor protein tyrosine kinases FGFR (1-4) (Qin et al., 2019). Each receptor contains an extracellular domain composed of three immunoglobulin (Ig)-like domains (domains I to III), a transmembrane domain, and an intracellular tyrosine kinase structure 2 (Turner and Grose, 2010). Fibroblast growth factors (FGFs) combined with FGFRs induce dimerization and phosphorylation of specific cytoplasmic tyrosine residues, and activate cytoplasmic signal transduction pathways (Carter et al., 2015). Studies have shown that genetic variation in the FGF/FGFR pathway is significantly related to the risk of ovarian cancer, treatment response and survival rate (Meng et al., 2014; Haugsten et al., 2010). Most studies have shown that the FGF/FGFR pathway plays a role in ovarian diseases, and it can be seen that the genes on this pathway are expressed in the ovaries. Erythrocyte membrane protein band 4.1 like 5 (EPB41L5) is a member of the band 4.1 superfamily (Hirano et al., 2008), locate on chromosome 2. It has a FERM domain at its N terminus and an actin-binding domain at its C terminus (Gamblin et al., 2018), which was related to cell movement. In addition, it can also participate in the development process of the body and can regulate the expression of downstream cadherin E (E-cadherin) in the process of epithelial cell-mesenchymal transition (EMT transformation) induced by transforming growth factor-β (TGF-β) (Li et al., 2019), the EMT transformation process played an important role in embryonic development, chronic inflammation, tissue reconstruction, cancer metastasis, and various fibrotic diseases. At the same time, EPB41L5 directly controlled the contractility and focal adhesion of actin, and regulated cell proliferation and migration (Hashimoto et al., 2016). The role of abnormal expression of EPB41L5 in malignant tumors has gradually been confirmed (Chen et al., 2019). It can be seen that the FGF/FGFR pathway has a certain influence in ovarian diseases. EPB41L5 was involved in embryonic development, so these two genes may be related to goat kidding performance. Based on the results of re-sequencing of Yunshang black goat (The Yunshang black goat is the first meat purpose black goat breed in China, which is bred from the Yunling goat as the dam and the Nubi goat as the sire (Tao et al., 2020a, b)) in the early stage of our laboratory, it is speculated that FGFR gene and EPB41L5 gene may have some influence on kidding performance. A MassARRAY® SNP typing technology was used to detect the genotypes of candidate genes and analyze the association between their polymorphisms and litter size traits. This study aimed to analyze the genetic polymorphisms of FGFR and EPB41L5 in Yunshang black goat, Jining grey goat and Liaoning cashmere goat population and their correlation with kidding performance, and provide more information about genetic markers for improving the prolificacy of goat.

A total of 768 goats (2~5 years old) were employed in this study, including 544 Yunshang black goats, 133 Jining grey goats, and 91 Liaoning cashmere goats. The higher prolificacy breeds included in this study was Jining Grey goats, and Liaoning Cashmere goats was regarded as lower prolificacy breeds (Wang et al., 2020; Chu et al., 2013). The breeding conditions of all experimental goats were the same. Yunshang Black goat had records of the number of litters born at least one litter, and some had records of litter weight at birth and litter weight at weaning (3 months old). Jugular blood samples from 768 goats of the three breeds were collected into EDTA-coated tubes and stored at –20°C until used for DNA isolation.

Samples of test goat were processed with the phenol–chloroform method for DNA isolation. The concentration of DNA samples were detected by using Nano Drop2000, and the quality of DNA samples were detected by using 1.5% agarose gel electrophoresis.

FGFR1 gene g. 12120297 A>G, EPB41L5 gene g. 64237881 A>C, g. 64266710 G>C and g. 64266715 G>C sites were used to detect by Sequenom MassARRAY® SNP typing technology. The related primer information was shown in Table I. The typing sample is DNA, the required amount of each sample is 20 μL, and the DNA concentration is 40~80 ng/μL.

A phylogenetic tree of EPB41L5 gene coding sequence in fifteen species (goat, pig, alpaca, human, horse, sheep, dog, cow, Norway rat, mouse, monkey, chicken, duck, fish and cat) were constructed with MEGA software. Sequences were obtained from the NCBI (https://www.ncbi.nlm.nih.gov/) Reference (NC_030809.1, NC_010457.5, NW_021964164.1, NC_000002.12, NC_009161.3, NC_056055.1, NC_051823.1, NC_037329.1, NC_051348.1, NC_000067.7, NC_041765.1, NC_052538.1, NC_051778.1, NC_007120.7, NC_018732.3, respectively). Prediction of the physical and chemical properties of EPB41L5 protein was performed using ProtParam. The secondary structure of EPB41L5 was analyzed using SOPMA. Prediction of the mRNA Structure of EPB41L5 was performed using RNAFOLD.

Based on genotyping data, genotype frequency, allele frequency, polymorphic information content (PIC), heterozygosity (He) and effective allele number (Ne) were calculated by using Microsoft Excel 2020 software. Hardy-Weinberg equilibrium test was carried out by Chi-square test. The one-way ANOVA was performed by IBM SPSS Statistics 22.0 software, and the LSD method was used for multiple comparisons. The correlation analysis between goat genotype and kidding phenotype data was carried out, and all data were expressed as mean±standard error.

The genotyping results showed that one SNP in FGFR1 (g.12120297 A>G) and three SNPs in EPB41L5 (g.64237881 A>C, g. 64266710 G>C and g. 64266715 G>C) were detected in all three breed of Yunshang Black goat, Jining Grey goat and Liaoning Cashmere goat. Two genotypes at FGFR1 gene g.12120297 A>G locus, AA and GA, were detected in 759 individuals with a detection rate of 98.8%. Three genotypes at EPB41L5 gene g. 64237881 A>C locus, AA, CC and CA, were detected in 756 individuals with a detection rate of 98.4%. The EPB41L5 g.64266710 C>G and g.64266715 G>C loci contained three genotypes, CC, GG and GC. The locus g.64266710 C>G was detected in 762 individuals with a detection rate of 99.2%, and the locus g.64266715 G>C in 748 individuals with a detection rate of 97.4% (Fig. 1).

It can be seen from Table II that the genotype frequency and gene frequency of EPB41L5 gene g.64237881 A>C, g.64266710 G>C and g.64266715 G>C locus differs between the high- and low-reproduction goat populations to a very significant level (P< 0.01). The dominant genotype of EPB41L5 gene g.64237881A>C locus in the high and low-reproduction populations was AA, and the dominant allele is A; the dominant genotype of EPB41L5 gene g. 64266710 G> C and g. 64266715 G> C locus in the high- and low-reproduction populations were GG, and the dominant allele was G. Population genetics statistics showed that there was a genetic polymorphism at FGFR1 gene g. 12120297 A> G locus in Yunshang Black goat, but there were no gene mutations in Jining Grey goat and Liaoning Cashmere goat. The FGFR1 gene g. 12120297A> G and EPB41L5 gene g. 64237881A> C were low polymorphisms in the three population (PIC< 0.25); EPB41L5 gene g. 64266710G> C and g. 64266715G> C were moderately polymorphic in Yunshang Black goat and Liaoning Cashmere goat (0.25< PIC< 0.5), and in Jining Grey goat were low polymorphism (PIC< 0.25). The χ2 test indicated that EPB41L5 gene g. 64237881A> C locus was under Hardy-Weinberg disequilibrium in Jining Grey goat (P< 0.05), and was under Hardy-Weinberg equilibrium in Yunshang Black goat and Liaoning Cashmere goat (P> 0.05). FGFR1 gene g. 12120297A> G, EPB41L5 gene g. 64266710G> C, g. 64266715G> C were under Hardy-Weinberg equilibrium in Yunshang Black goat, Jining Grey goat and Liaoning Cashmere goat (P> 0.05).

Association analysis indicated that there were no significant correlations between FGFR1 gene g.12120297A>G locus and EPB41L5 gene g.64266715G>C locus and litter size, litter weight at birth and weaning (P> 0.05). There was significant correlation between EPB41L5 gene g.64237881A>C and g.64266710 G>C and litter size (P< 0.05), and there was no significant correlation between the two loci and litter weight at birth (P> 0.05). There was no significant correlation between EPB41L5 gene g.64237881A>C and litter weight at weaning (P> 0.05), there were significant correlation between EPB41L5 gene g.64266710C>G and litter weight at weaning (P< 0.05).

A phylogenetic tree of EBP41L5 gene coding sequence for other animal species was constructed. The results showed that there was the highest homology between goats and sheep, and there was the higher similarity between goats and cattle, and the sequences of the other species showed relatively large differences between them (Fig. 2).

Note: F, Upstream primer; R, Downstream primer; E, Extension primer.

Note: P< 0.05 means the difference is significant; P< 0.01 means the difference is extremely significant.

It can be seen from Table V that the molecular weight, formula, total number of atoms, instability index, aliphatic index and grand average of hydropathicity of the protein before and after the mutation have happened (Table V). It can be seen from the figure that the hydrophilic amino acids are more than the hydrophobic amino acids before and after EPB41L5 mutation (Fig. 3A, 3B). Therefore, the protein was hydrophobic before and after the mutation, but grand average of hydropathicity has changed. It can be seen that the physicochemical properties of the protein are changed due to the mutation of the three sites, which may cause the function of the protein to change, thereby affecting the kidding traits of the Yunshang Black goat.

Note: PIC, HE and NE represent polymorphic information content, heterozygosity and effective number of alleles, respectively; numbers in the parentheses represent the number of detected sheep of each genotype; YS, JN and LN represent Yunshang Black goat, Jining Grey goat, Liaoning Cashmere goat, respectively; P> 0.05 indicates the locus was under Hardy-Weinberg equilibrium; P< 0.05 indicates the locus was under Hardy-Weinberg disequilibrium.

Note: Numbers in the parentheses next to litter size represent the number of sheep of each genotype; Different lower-case letters in the same group indicate significant difference (P < 0.05).

For two missense mutations found in EPB41L5 of goat, the secondary structure of the protein and mRNA before and after the mutation at g. 64266710 G> C, g.64266715 G> C and g. 64237881 A>C were predicted, respectively. It can be seen that alpha helix, extended strand, beta turn and random coil in the secondary structure of the protein have changed before and after the mutation (Table VI). The results showed that the secondary structure of the protein (Fig. 4A, B) and mRNA (Fig. 4C, D) did change significantly after the mutation. And its minimum structural free energy is reduced, from -716.13kJ mol-1 to -719.88.kJmol-1.

Fibroblast growth factor receptor 1 (FGFR1) was located on goat chromosome 27. The binding of FGF led to FGFR dimerization and conformational changes, which phosphorylate the tyrosine kinase domain. Fibroblast growth factor receptor substrate 2 was a crucial adaptor protein in the FGFR signaling cascade. The phosphorylated tyrosine residues in turn became its docking site, and its phosphorylation led to adaptor proteins such as the recruitment of SOS and GRB2, activated downstream MAPK, PI3K/AKT and other signaling pathways (Datta et al., 2017). The physiological roles of FGF/FGFR signaling pathway included embryonic development, wound healing, angiogenesis and tissue crosstalk (Heinzle et al., 2011). In addition, the FGF/FGFR signaling pathway was essential for cell proliferation, migration, differentiation, apoptosis and survival (Dieci et al., 2013). Studies have shown that abnormal FGF/FGFR signals are closely related to the pathogenesis of a variety of cancers, including lung cancer, gastric cancer, breast cancer, ovarian cancer, urothelial cancer and endometrial cancer (Dorey and Amaya, 2020). Erythrocyte membrane protein band 4.1 like 5 (EPB41L5) also known as BE37, YMO1, had a molecular weight of 85 kD. EPB41L5 discovered for the first time in November 2007, which was a protein related to cytoskeleton and cell movement (Rhim et al., 2012), and it was widely expressed in mammalian embryonic mesoderm and endoderm epithelial tissues (Tepass, 2009). It had a FERM domain at their N-terminus, but lacks an actin binding site at their C-terminus (Hashimoto et al., 2016). EPB41L5 was involved in the regulation of many important biological processes, included the cellular actin skeleton, cell adhesion molecules, EMT, basement membrane formation, cell polarity formation and other processes. Research has shown that the FERM domain contained in EPB41L5 can cross-link with the myosin protein containing the myosin tail homology 4 (MyTH4) domains. Involved in regulating the formation of pseudopodia and cell movement, and can be combined with the cytoskeleton protein Talinl, PTK, etc., (Sordella, 2003). Hirano’s (Hirano et al., 2008) research showed that EPB41L5 overexpression can reduce Ecadherin on the cell membrane surface through the action of its FERM domain, thereby promoting the EMT process. Gosens et al. (2007) found that CRB-MPP5-EPB41L5 complex can regulate a variety of cell adhesion factors, it played an important role in the tight junctions between cells and in the formation of cell polarity, and was closely related to epithelial cell polarity and basement membrane formation. In addition, experiments have shown that EPB41L5 can affect the expression of cell adhesion molecules and skeletal proteins in embryonic endoderm epithelial cells (Hirano et al., 2008). At the same time, EPB41L5 can regulate cell polarity and was an essential protein for basement membrane formation and embryonic gastrulation (Takakuwa et al., 1986). In summary, the expression of EPB41L5 is related to embryonic development and may affect goat kidding traits. The genotype frequency and allele frequency of EPB41L5 reached extremely significant levels between the high- and low-reproduction goat populations (P< 0.01). Therefore, it is speculated that its mutation may have a certain impact on fecundity, and it needs to be further verified by experiments. Association analysis showed that there were significant correlation between EPB41L5 gene g.64237881A> C and g.64266710G> C and number of litter size (P< 0.05). There were significant correlation between EPB41L5 gene g. 64266710C> G and litter weight at weaning (P< 0.05). It indicates that EPB41L5 gene g. 64237881A> C and g. 64266710G> C may be critical sites for the number of litter size and litter weight at weaning. Therefore, in order to provide a reference for goat molecular breeding, this study aimed to analyze the genetic polymorphisms of FGFR and EPB41L5 in Yunshang Black goat, Jining Grey goat and Liaoning Cashmere goat population and their association with kidding performance.

Correlation analysis showed that g. 64237881A> C, g. 64266710G> C of EPB41L5 gene may affect kidding performance. In order to further explore the changes in the structure of the protein before and after the mutation of g. 64237881A> C and g. 64266710G> C, two sites were analyzed by bioinformatics. Bioinformatics analysis of EPB41L5 gene showed that the gene conforms to the laws of biological evolution. Before and after the mutation, the physical and chemical properties, the secondary structure of the protein and mRNA were changed. Therefore, it is possible that mutations in the site may result in changes in the structure of the protein, which in turn may result in changes in function. Therefore, g. 64237881A> C, g. 64266710G> C of EPB41L5 gene may be the key site for litter size and litter weight at weaning.

There was significant correlation between EPB41L5 gene g. 64237881A> C and g. 64266710G> C and litter size (P< 0.05). There was significant correlation between EPB41L5 gene g. 64266710C> G and litter weight at weaning (P< 0.05). Therefore, these results suggested that the EPB41L5 gene g. 64237881A> C and g. 64266710G> C loci are suitable as a molecular marker for litter size in Yunshang Black goat, and g. 64266710G> C locus is suitable as a selection marker for litter weight at weaning.

This research was financially supported by China Agriculture Research System of MOF and MARA (CARS-38), Major Science and Technology Project of Yunnan Province (202102AE090039), the Basic Research Foundation Key Project of Yunnan Province (202001AS070002), Agricultural Science and Technology Innovation Program of China (CAAS-ZDRW202106 and ASTIP-IAS13), Key Research and Development Project of Shanxi Province (201903D211008), Shanxi Province 1331 Project Key Disciplines of Animal Sciences (J201911301), Shanxi Agricultural University Doctoral Research Initiation Project (2021BQ04).

All authors have declared no conflicts of interest.