Polymorphisms of the Growth Hormone Releasing Hormone Receptor Gene Affect Body Conformation Traits in Chinese Dabieshan Cattle

Simple Summary Increasing growth traits can impact the economic benefit of livestock. Polymorphisms can offer us potential molecular markers for the growth traits of cattle, and provide a reference for cultivating an improved cattle population according to the markers in the future. In this study, the growth hormone-releasing hormone receptor gene polymorphism is significantly associated with the body conformation traits of Chinese Dabieshan cattle, and it could be used as a candidate genetic marker for body conformation trait selection in cattle. This study further verified that the GHRHR gene plays a vital role in Chinese cattle, and provides a data basis for the research and improvement of cattle breeding programs in China. Abstract This study was performed to expose the polymorphisms of the growth hormone-releasing hormone receptor gene in Chinese Dabieshan cattle, evaluate its effect on body conformation traits, and find potential molecular markers in Chinese cattle. The GHRHR structure and the phylogenetic tree were analyzed using bioinformatics software. The polymorphism of the GHRHR gene in 486 female cattle was genotyped by PCR-RFLP and DNA sequencing, and the association between SNPs and body conformation traits of Chinese Dabieshan cattle was analyzed by one-way ANOVA in SPSS software. GHRHR was often conserved in nine species, and its sequence of cattle was closest to sheep and goats. Six polymorphic SNPs were identified, g.10667A > C and g.10670A > C were missense mutation. The association analysis indicated that the six SNPs significantly influenced the body conformation traits of Chinese Dabieshan cattle (p < 0.05). Six haplotypes were identified and Hap1 (-CAACGA-) had the highest frequency (36.10%). The Hap3/5 (-GCCCCCGGAAGG-) exhibited a significantly greater wither height (WH), hip height (HH), heart girth (HG), and hip width (HW) (p < 0.05). Overall, the polymorphisms of GHRHR affected the body conformation traits of Chinese Dabieshan cattle, and the GHRHR gene could be used as a molecular marker in Dabieshan cattle breeding programs.


Introduction
Dabieshan cattle area Chinese native breed and are widely rearedin the middle and lower reaches of the Yangtze River [1]. Based on the mitochondrial DNA analysis and Y-SNPs and Y-STRs markers tests, Dabieshan cattle belong to the Southern cattle type in China, have two types of maternal origin from Bos taurus and Bos indicus, and the displacement loop regions had rich genetic diversity (π = 0.023) [2,3]. Compared with exotic commercial cattle breeds, Dabieshan cattle have some important traits, including adaptation to low-quality feed resources, fine fat deposition capabilities, wet-heat tolerance and disease resistance, whereas slow growth was prioritized for genetic improvement. The

Ethics Statement
The experimental animals and procedures performed in this study were approved by the Animal Care and Use Committee of the Anhui Academy of Agricultural Sciences (approval number A19-CS20).

Phenotypic Data and DNA Sample Collection
A total of 486 female Dabieshan cattle aged 24 to 30 months were randomly chosen from National Species Resources Protection Farm (Anqing, China). According to NRC Animals 2022, 12, 1601 3 of 13 standards (Nutrient Requirement of Beef Cattle), the animals on the farm were fed a total mixed ration (TMR), which comprised 25% concentrate and 75% roughages of dry straw and corn silage, and water was offered ad libitum.
The phenotype performance for the body conformation traits in Dabieshan cattle, including body length (BL), wither height (WH), hip height (HH), heart girth (HG), abdominal girth (AGR), hip width (HW), and pin bone width (PBW), were collected following the methods of Wang et al. [18] and Yang et al. [19]. Ear marginal tissues collected from these animals were used to extract genomic DNA according to the TIANamp Genomic DNA Kit (TIANGEN, Beijing, China) procedure. Genomic DNA were measured with the spectrophotometer, and then diluted to 50 ng/µL for polymerase chain reaction (PCR) analysis.

Primers, Polymerase Chain Reaction, and Gel Electrophoresis
Primer Premier 5 software (PREMIER Biosoft International, San Francisco, CA, USA) was used to design the primers (Table 1) for amplifying the bovine GHRHR gene (GenBank accession number: NC_037331.1). The fragments were amplified with 25 µL of reaction volume, which contained 1.0 µL of each forward and reverses primers (100 ng/µL), 1 µL genomic DNA (50 ng/µL), 12.5 uL of 2 × Taq Mix (TIANGEN, Beijing, China), and 9.5 µL double distilled water (ddH 2 O). The PCR amplification was performed according to the procedure, including pre-denaturation at 94 • C for 5 min, 32 cycles of denaturation at 94 • C for 30 s, annealing at 60 • C for 30 s, and extension at 72 • C for 30 s, followed by a final extension at 72 • C for 5 min. The PCR products of F1R1 were genotyped by PCR-RFLP. The digestion reaction mixture contained 5 uL PCR products, 1 uL 10 × buffer, 10 U restriction enzyme (PshAI, NEB, Ipswich, MA, USA), and sterile water. Samples were incubated at 37 • C for restriction enzyme overnight, and the digested products were analyzed by 1.5% agarose gel electrophoresis, stained with GeneGreen Nucleic Acid Dye (TIANGEN, Beijing, China) in 1% TBE. The PCR products of F2R2 were also electrophoresed on a 1.5% agarose gel, and then purified and sequenced through Sangon Biotech (Shanghai, China). GHRHR-F for forward and GHRHR-R for reverse.

Statistical Analysis
The allelic and genotypic frequencies of six SNPs were measured by direct counting.Hardy-Weinberg equilibrium (HWE) was assessed using POPGENE 3.2 software package with the Chi-squared (χ 2 ) test. Population genetic indexes, such as gene heterozygosity (H e ), gene homozygosity (H o ), effective allele numbers (N e ) and polymorphism information content (PIC), were analyzed according to Chakraborty and Nei [20]. LD plot was performed to analyze the Linkage disequilibrium (LD)of SNPs and haplotype block of Haploview software was used to obtain the haplotypes [21,22]. The diplotype was assessed according to the SNPs and haplotypes.
The association analysis was calculated by one-way ANOVA in the SPSS software (version 24.0) (Armonk, NY, USA), which was used in our previous study [23]. The general linear model was used to determine the association between SNPs and body measurement traits. The equation was shown as follows: y ijk = u + g i + a j + s k + e ijk (1) where y ijk was the phenotypic observation, u was the mean, g i was the fixed effect of genotype, a j was the random effect of sire, s k was the fixed effect of age, and e ijk was the residual effect. The results were presented as the mean ± standard error, and p < 0.05 was considered a significant difference.

Species Homology, Phylogenetic Tree, and Conservative Estimation
The cDNA of the GHRHR in cattle consisted of 102-bp 5 UTR sequences, followed by a 1326-bp open reading frame and 993-bp 3 UTR. A total of 441 amino acid residues were encoded in the coding region with a molecular weight of 49.25 kD and an isoelectric point of 6.37. In Table 2, both the 5 -donor and 3 -acceptor splice sites in the GHRHR conformed to the GT-AG rule. Figure 1 shows the conservative analysis of the GHRHR protein in nine animals (human, rat, mouse, cattle, pig, goat, sheep, horse, and chicken), and signal peptide, HormR, and Pfam were often observed. The neighbor-joining phylogenetic tree ( Figure 2) revealed that cattle were closest to goats and sheep for the GHRHR sequence, while humans, rats, mice, and chickens were fartherfrom the cattle branch. In total, 12 significant motifs were found using MEME online in the super secondary protein structure of the GHRHR gene ( Figure 3), indicating that it plays a similar role in the nine animals. eral linear model was used to determine the association between SNPs and body measurement traits. The equation was shown as follows: yijk = u + gi + aj + sk + eijk (1) Where yijk was the phenotypic observation, u was the mean, gi was the fixed effect of genotype, aj was the random effect of sire, sk was the fixed effect of age, and eijk was the residual effect. The results were presented as the mean ± standard error, and p< 0.05 was considered a significant difference.

Gene Number Exon Size(bp) Intron Size (bp) 5′Splice Donor 3′Splice Acceptor
The neighbour-joining phylogenetic tree of GHRHR protein sequences in nine animals. Figure 2. The neighbour-joining phylogenetic tree of GHRHR protein sequences in nine animals.

Identification of SNPs in GHRHR
In this study, six SNPs including g.2052C > G in intron 1, g.10667A > C and g.10670A > C in exon 10, and g.10684C > G, g.10699A > G, and g.10763A > G in intron 10, were identified in the GHRHR by PCR-RFLP and sequencing. The agarose gel electrophoresis of PCRamplified products using F1R1 and F2R2 primers were shown in Figure 4a,b, and the PCR-PshAI-RFLP analysis at g.2052C > G locus was present in Figure 4c. The 498 bp PCR fragments of F2R2 primers can be digested by PshAI into 289 bp and 209 bp for allele C. The sequenced map of five SNPs in the GHRHR was shown in Figure 5. All these SNPs are biallelic.
12, x FOR PEER REVIEW 6 of 13 Figure 3. The significant motifs of GHRHR protein sequences in nine animals. These colors were given by MEME suit system during analysis.

Identification of SNPs in GHRHR
In this study, six SNPs including g.2052C > G in intron 1, g.10667A > C and g.10670A > C in exon 10, and g.10684C > G, g.10699A > G, and g.10763A > G in intron 10, were identified in the GHRHR by PCR-RFLP and sequencing. The agarose gel electrophoresis of PCR-amplified products using F1R1 and F2R2 primers were shown in Figure 4a,b, and the PCR-PshAI-RFLP analysis at g.2052C > G locus was present in Figure 4c. The 498 bp   The SNPs in exon10 were missense mutations (g.10667A > C: Glutamine (Gln)→ Proline (Pro); g.10670A > C: Histidine (His)→Pro). Figure 6 shows the predicted 3D structure and the value of quality model energy analysis (QMEAN) of the normal and mutant GHRHR (g.10667A > C and g.10670A > C). The figures indicated there was no change in   The SNPs in exon10 were missense mutations (g.10667A > C: Glutamine (Gln)→ Proline (Pro); g.10670A > C: Histidine (His)→Pro). Figure 6 shows the predicted 3D structure and the value of quality model energy analysis (QMEAN) of the normal and mutant GHRHR (g.10667A > C and g.10670A > C). The figures indicated there was no change in The SNPs in exon10 were missense mutations (g.10667A > C: Glutamine (Gln)→ Proline (Pro); g.10670A > C: Histidine (His)→Pro). Figure 6 shows the predicted 3D structure and the value of quality model energy analysis (QMEAN) of the normal and mutant GHRHR (g.10667A > C and g.10670A > C). The figures indicated there was no change in the 3D structure (Figure 6a-d), whereas the QMEAN value of the mutant-type was significantly different from the wild-type (Figure 6e-h). As shown in Figure 6e,f, the mutation of Gln to Pro at the 323 sites led to the value of QMEAN increasing from −4.70 to −4.73, however, the torsion and all atoms of the wild-type had a lower level than the mutant-type, while the solvation and Cβ were opposite. Figure 6e,g show that the conversion of His to Pro at the 324 sites leads to a decrease in the value of QMEAN, solvation, torsion, Cβ, and all atoms. Additionally, the conversion of both Gln and His to Pro caused the rise of all atoms and torsion, and reduction of QMEAN, Cβ and solvation (Figure 6e,h). the 3D structure (Figure 6a-d), whereas the QMEAN value of the mutant-type was significantly different from the wild-type (Figure 6e-h). As shown in Figure 6e,f, the mutation of Gln to Pro at the 323 sites led to the value of QMEAN increasing from −4.70 to −4.73, however, the torsion and all atoms of the wild-type had a lower level than the mutanttype, while the solvation and Cβ were opposite. Figure 6e,g show that the conversion of His to Pro at the 324 sites leads to a decrease in the value of QMEAN, solvation, torsion, Cβ, and all atoms. Additionally, the conversion of both Gln and His to Pro caused the rise of all atoms and torsion, and reduction of QMEAN, Cβ and solvation (Figure 6e, h).

Genetic Diversity, and Hardy Weinberg Equilibrium
In the Dabieshan cattle, 486 female cattle were genotyped by PCR-RFLP and sequencing. The genetic diversity including genotypic, allelic frequency and genetic indexes of the six detected SNPs in the GHRHR was evaluated in Table 3. The G allele at g.10699A > G locus was most prevalent (0.71), whereas the GG genotype at g.10684C > G was more frequent than others. The values of heterozygosity (He) were 0.414~0.499 and effective allele numbers (Ne) values were1.707~2.000. Polymorphism information content (PIC) values ranged from 0.328 to 0.375, indicating that the genetic diversity of Dabieshan cattle at these six SNP sites was at a medium level (0.25 < PIC < 0.5).

Genetic Diversity, and Hardy Weinberg Equilibrium
In the Dabieshan cattle, 486 female cattle were genotyped by PCR-RFLP and sequencing. The genetic diversity including genotypic, allelic frequency and genetic indexes of the six detected SNPs in the GHRHR was evaluated in Table 3. The G allele at g.10699A > G locus was most prevalent (0.71), whereas the GG genotype at g.10684C > G was more frequent than others. The values of heterozygosity (H e ) were 0.414~0.499 and effective allele numbers (N e ) values were1.707~2.000. Polymorphism information content (PIC) values ranged from 0.328 to 0.375, indicating that the genetic diversity of Dabieshan cattle at these six SNP sites was at a medium level (0.25 < PIC < 0.5).

Linkage Disequilibrium, and Haplotype Analysis
Chi-square tests (p < 0.05) were calculated for the Hardy-Weinberg equilibrium (HWE) in detected SNPs. Table 3 shows that six SNPs except forg.10763A > G were in the Hardy-Weinberg equilibrium. Table 4 presents the linkage disequilibrium among the detected SNPs, the D' values were in the range of 0.216 and 1.000, and r 2 values were on a scale of 0.012 to 0.919. Several studies indicated that r 2 is more suitable for measuring LD because of its insensitivity to allele frequencies [24]. The r 2 values in this study showed that g.10670A > C, g.10684C > G, and g.10763A > G were in a relatively strong linkage, and g.10667A > C was relatively strongly linked with g.10699A > G in Dabieshan cattle (r 2 > 0.33). Table 4. Linkage equilibrium analysis between SNP markers in GHRHR. The haplotype analysis at the six loci of the GHRHR gene was illustrated in Table 5. Six dominating haplotypes were identified in the Dabieshan cattle, Hap1 (-CAACGA-) was the most frequent (36.10%), followed by Hap2 (-GACGGG-), Hap3 (-CCCGAG-), and Hap4 (-GAACGA-). Hap5 (-GCCGAG-) and Hap6 (-CACGGG-) showed the lowest frequency in our analysis. Table 5. Haplotypes analysis of the six SNPs in GHRHR.

Association Analysis
The association between diplotype and body conformation traits was analyzed in Table 7. In Dabieshan cattle, animals carrying the Hap3/5 (-GCCCCCGGAAGG-) exhibited a significantly greater WH, HH, HG, and HW than the others (p < 0.05). The frequencies of diplotypes < 5.0% were not considered. which may be due to the population size and artificial selection. The r 2 -values, among the g.10670A > C, g.10684C > G, and g.10763A > G sites, and between g.10667A > C and g.10699A > G sites, were sufficiently strong for mapping. The results may be caused by the lower recombination and higher genotypic variation at these sites [29].
There were many agreeing reports regarding the effect of the GHRHR mutations on the body conformation traits of cattle in previous studies. Zhang et al. [14] indicated that a mutation in 5 UTR of the GHRHR had a significant effect on the Nanyang cattle phenotype. In this study, there was a significant association between the six identified SNPs and body conformation traits in Dabieshan cattle. Moreover, Hap3/5 (-GCCCCCGGAAGG-) had a greater WH, HH, HG, and HW, the Hap3/3 (-CCCCCCGGAAGG-) and Hap2/5 (-GGACCCGGAGGG-) exhibited greater HG and HW than animals with other diplotypes, which indicated that Hap3 (-CCCGAG-) and Hap5 (-GCCGAG-) might result in an increasing body conformation trait in Chinese Dabieshan cattle.
Interestingly, the SNPs g.10667A > C and g.10670A > C in exon 10 were missense mutations (Gln → Pro; His → Pro). Figure 6 indicates that the QMEAN value, all atoms, and torsion of mutant 1 were higher than the wild-type and other mutants, indicating that the protein of mutant 1 might have a stronger potential energy and torsion angle in connection with pairwise atom distances than other types. The Cβ and solvation of mutant 3 were lower than those of other types, which suggest that the mutant 3 protein might exhibit a lower residue embedding ability and a weaker C-β potential energy of interaction. In accordance with the study of Wang et al. [30], the SNPs located in intron 1 and 10, show that the non-coding region variation in the GHRHR significantly affects the growth and development of Dabieshan cattle. Increasingly, the literature has shown that intron mutation might affect the correct processivity of the mRNA by disrupting the splice site, altering the secondary structure, affecting protein function to regulate the gene expression level or determine tissue-specific expression patterns [31].
In general, six SNPs including g.2052C > G, g.10667A > C, g.10670A > C, g.10684C > G, g.10699A > G, and g.10763A > G of the GHRHR gene were genotyped in Chinese Dabieshan cattle. The six detected SNPs played a significant role in the body conformation traits, and animals with Hap3/5 exhibited significantly improving body conformation traits, meaning that the GHRHR could be used as a marker in Dabieshan cattle breeding programs.

Conclusions
This study investigated six novel SNPs in the GHRHR gene in Chinese Dabieshan cattle. Our results show that all six loci contributed significantly to the body conformation traits of cattle (p < 0.05), and Hap3/5 was attributed to desirable body characteristics. These findings indicated that the GHRHR could be used as a DNA marker for improving the performance of Dabieshan cattle.  Institutional Review Board Statement: The animal study protocol was approved by the Animal Care and Use Committee of the Anhui Academy of Agricultural Sciences (approval number A19-CS20 and date of approval 16 February 2020).

Informed Consent Statement: Not applicable.
Data Availability Statement: All data is provided in the manuscript.