Polymorphism detection of PRKG2 gene and its association with the number of thoracolumbar vertebrae and carcass traits in Dezhou donkey

Background Previous studies have shown that the protein kinase cGMP-dependent 2 (PRKG2) gene is associated with dwarfism in humans, dogo Argentines, and Angus cattle, as well as with height and osteoblastogenesis in humans. Therefore, the PRKG2 gene was used as the target gene to explore whether this gene is associated with several thoracolumbar vertebrae and carcass traits in Dezhou donkeys. Results In this study, fifteen SNPs were identified by targeted sequencing, all of which were located in introns of the PRKG2 gene. Association analysis illustrated that the g.162153251 G > A, g.162156524 C > T, g.162158453 C > T and, g.162163775 T > G were significantly different from carcass weight. g.162166224 G > A, g.162166654 T > A, g.162167165 C > A, g.162167314 A > C and, g.162172653 G > C were significantly associated with the number of thoracic vertebrae. g.162140112 A > G was significantly associated with the number and the length of lumbar vertebrae, and g.162163775 T > G was significantly associated with the total number of thoracolumbar vertebrae. Conclusion Overall, the results of this study suggest that PRKG2 gene polymorphism can be used as a molecular marker to breed high-quality Dezhou donkeys. Supplementary Information The online version contains supplementary material available at 10.1186/s12863-022-01101-6.


Introduction
The donkey industry is an integral part of modern animal husbandry, significantly increasing the economic income of both free-range farmers and large farms. Donkey meat is delicious food consumed in some countries, and is highly nutritious and has a unique flavor [1]. Donkeys are uniparous animals and have long growth cycles. Dezhou donkeys reach sexual maturity at about 12-15 months, so molecular breeding of donkeys to improve meat production is necessary and urgent. The number of thoracic vertebrae ranged from 17 to 19, and the number of lumbar vertebrae ranged from 5 to 6 in Dezhou donkey [2]. Previous studies have found that changes in the number of thoracolumbar vertebrae can provide economic benefits. An extra vertebra increases carcass weight by 6 kg [2]. Therefore, it is of great significance to breed multiple thoracolumbar donkeys to improve the quantity of meat.
Many studies have previously demonstrated that variation in the number of thoracolumbar vertebrae can lead to changes in economic traits such as body length and carcass weight in pigs [3] and sheep [4]. In recent years, selection and breeding for multiple thoracolumbar vertebrae traits in pigs, cattle, and sheep have been carried out to analyze the primary loci for thoracolumbar numbers. A point mutation in intron 4 of the ActRIIB gene in Small Tailed Han sheep was associated with variation in vertebral number [5]. The TGFβ3 gene was a candidate gene for the number of vertebrae traits in pigs. The g.105179474 G > A mutation locus on chromosome 7 was associated with the number of ribs and thoracolumbar vertebrae [6]. g.19034 A > C locus of VRTN gene can be used as a potential molecular marker for multiple thoracic vertebrae number in Beijing black pigs [7]. However, the selection and breeding for multiple thoracolumbar vertebrae in donkeys have just started. In donkey, the HOXC8 g.15179224C > T was significantly associated with lumbar vertebrae length (P < 0.05), and the g.15179674G > A locus was shown to be significantly associated with the number of lumbar vertebrae (P < 0.05) [8]. The NR6A1 g.18114954C > T is significantly associated with lumber vertebrae number and the total number of thoracolumbar, and individuals with TT genotype had significantly larger value than CC genotype (P < 0.05) [2]. Therefore, it is valuable and essential to identify genes affecting multiple thoracolumbar vertebrae numbers and carcass traits in Dezhou donkeys.
Many studies have shown that the PRKG2 gene was associated with growth traits and skeletal development. PRKG2 gene is located on chromosome 3 in donkeys and contains eighteen exons and seventeen introns [9] (Fig. 1). Studies have demonstrated that the PRKG2 gene was associated with dwarfism in American Angus cattle [10], dogo Argentines [11], and humans [12]. The PRKG2 gene was identified as a candidate gene for human height by genome-wide analysis for copy number variants (CNVs) of 162 patients (149 families) with short stature [13]. The previous studies indicated the PRKG2 gene as a candidate for osteoblastogenesis [14]. Considering that the PRKG2 gene affects human height, human height is equivalent to donkey body length, and donkey body length is related to the number of thoracolumbar vertebrae, so the PRKG2 gene assumes to be associated with the number of thoracolumbar vertebrae and carcass traits. However, the association of the PRKG2 gene with the number of thoracolumbar vertebrae and carcass traits in Dezhou donkeys has not been reported.
In the present research, genetic variation in the PRKG2 gene of the Dezhou donkey has been studied using targeted sequencing technology. The targeted sequencing method is a technique to achieve accurate genotype detection by high-depth resequencing of target genes, which has the advantages of high stability, tolerance to sequence conservation and GC content, and can achieve excellent capture efficiency with flexible marker types and capture types [15]. Currently, targeted sequencing technology is widely used in human [16], plant [17], and animal [18]. The study aimed to investigate the genetic variation of the PRKG2 gene and its correlation with number of thoracolumbar vertebrae and carcass traits in Dezhou donkeys, and provide a specific theoretical basis for molecular breeding of Dezhou donkeys.

Ethics statement
The experimental animals and methods used in this study were approved by the Animal Policy and Welfare Committee of Liaocheng University (No. LC2019-1). The care and use of laboratory animals fully comply with local animal welfare laws, guidelines and policies.

Animals and phenotypes
Blood samples and trait data were collected from 406 2-year-old Dezhou donkeys at a slaughterhouse in Dezhou, Shandong Province. The 406 Dezhou donkeys in this study were all males and had the same feeding environment. Blood samples were collected from the jugular vein of donkeys using EDTA blood collection tubes and stored in a -20 °C refrigerator immediately. The relevant body size traits of donkeys were measured and recorded. Body height, body length, and chest circumference were measured under the National Standard of the People's Republic of China, "Dezhou Donkey." Carcass weight, the number of lumbar vertebrae, the number of thoracic vertebrae, the length of lumbar vertebrae, the length of thoracic vertebrae, the total number of thoracic and lumbar vertebrae were measured after humanely slaughtered. Carcass traits and the number of thoracolumbar vertebrae data were collected according to the method of Liu et al. (2022). All measurements are performed by the same operator to reduce human error. Table S1 is a summarizes the number of thoracolumbar vertebrae and carcass traits of 406 donkeys. Table S2 shows the mean overall situation of donkeys' thoracolumbar number and carcass traits, and the value is Means ± SE.

DNA extraction
Genomic DNA was extracted from blood samples using the TIANamp Blood DNA Kit (Tiangen, Beijing, China). After extraction, genomic DNA concentration was measured using a spectrophotometer (B500, Metash, China); a working solution was prepared and adjusted to 30 ng/ µL. The samples were placed in a − 20 °C refrigerator for later use.

SNP detection and genotyping
The 406 genomic DNA samples were sent to Molbreeding Biotechnology Co., Ltd. (Shijiangzhuang, China) for genotyping of the PRKG2 gene by Targeted Sequencing.
A total of 1292 probes were used in the targeted sequencing, covering 92.39% of the PRKG2 gene with reference sequence of the donkey PRKG2 gene (assembly ASM1607732v2; NC_052179; GCA_016077325.2). SNPs with genotype frequencies less than 5% in targeted sequencing results were removed.

SNPs validation
Sanger sequencing was used to verify the results of targeted sequencing. SNPs located at genomic position 162,150,000-162,160,000 bp in chromosome 3 were randomly selected for validation by Sanger sequencing, and the mutation sites in this region included g.162153251 G > A, g.162156524 C > T, and g.162158453 C > T. Three pairs of primers were designed to amplify three selected SNPs (g.162153251 G > A, g.162156524 C > T, g.162158453 C > T) in the PRKG2 gene using Primer Premier 5.0 software ( Table 1). The PCR amplification was performed in a total of 25 μL reaction, 12. The cycling parameters were as follows: pre-denaturation at 96 ℃ for 5 min, denaturation at 96 ℃ for 20 s, annealing at 62 ℃ for 30 s, and extension at 72 ℃ for 30 s. Each subsequent cycle is reduced by 1 ℃ until 52 ℃, for 10 cycles. 20 s of denaturation at 96 ℃, 30 s of annealing at 52 ℃, and 30 s of stretching at 72 ℃, 35 cycles. 10 min of extension at 72 ℃. 4 ℃ of storage. The specificity of the PCR products was detected using a 2% agarose gel, and samples that were detected for specificity and correct product size were sent to BGI Genomics Co., Ltd (Shanghai, China) for Sanger sequencing, and the results were analyzed using Chromas software (Version V2.6.5, Technelysium Pty Ltd., Queensland, Australia).

Statistical analyses
Genotype frequencies, allelic frequencies, and the Hardy-Weinberg equilibrium (HWE) were examined using Excel. Population genetic parameters, including homozygosity (Ho), heterozygosity (He), effective allele number (Ne) and the polymorphism information content (PIC) were analyzed using online software (http:// www. msrca ll. com/, accessed on 24 March 2022) [19]. The association of fifteen SNPs and haplotype combinations of the PRKG2 gene with the thoracolumbar number and carcass traits was analyzed using a general linear model where Y is the individual phenotypic measurements, µ represents the mean for each trait, a represents the fixed factor genotype, e represents the random error. Least squares means with standard errors were used for the different genotypes and for the number of thoracolumbar vertebrae as well as the carcass traits. Multiple comparisons of the associations were based on Bonferroni-corrected p-Values. The different genotypes were considered as fixed effects, the random error as a random effect and the number of thoracolumbar vertebrae and carcass traits as the dependent variable [21]. Linkage disequilibrium (LD) and haplotype construction were performed using Haploview 4.1 [22], and haplotypes with frequencies greater than 0.05 were constructed.

SNPs identification and genotyping
Targeted sequencing results showed that a total of 485 SNPs were identified (Table S3). Among them, 11 SNPs were located in exons, 457 SNPs were located in introns, 17 SNPs were located downstream of PRKG2 gene. However, 470 SNPs had a genotype frequency of less than 5%, Y ij = µ + ai + eij therefore statistics will not been applied to these data. The locations of these fifteen SNPs are shown schematically in Fig. 1. These fifteen SNPs of PRKG2 gene were genotyped using sequencing, which generated three genotypes for all locus. The genotyping results of fifteen SNPs of PRKG2 gene are shown in Table S4. The Sanger sequencing results of the three SNPs (g.162153251 G > A, g.162156524 C > T and g.162158453 C > T) were consistent with the targeted sequencing results. Three samples were randomly selected at three sites from 406 Dezhou donkey DNA samples were randomly selected as the amplification template for three SNPs, and the amplification products were added into 1% agarose gel for electrophoresis identification. Electrophoresis results showed that the bands were single, clear and bright, in line with the expected fragment size.

Genetic parameter analysis
The genotype and allele frequency were calculated ( Table 2). The mutant allele frequency of g. 162,140,112 A > G was the highest, and the normal allele frequency of g. 162,153,251 G > A was the highest. g.162153251 G > A, g.162156524 C > T and g.162216538 G > A were not in HWE. The values of Ho for the fifteen SNPs ranged from 0.2705 to 0.7333, He for the fifteen SNPs ranged from 0.2667 to 0.7295, and Ne for the fifteen SNPs ranged from 1.3636 to 3.6966. Only g.162153251 G > A was in low polymorphism (PIC < 0.25), while the other mutation sites were in moderate polymorphism

Linkage disequilibrium analysis and haplotype construction
Linkage disequilibrium (LD) analysis of the remaining loci showed a strong association between every two SNPs (r 2 > 0.33) (Fig. 2). Block 1 consisted of two SNPs (g.162158453 C > T, g.162160146 T > C). In block 1, the linkage disequilibrium of g.162158453 C > T with g.162160146 T > C was not very strong (r 2 < 0.33). Block 2 consisted of six SNPs (g.162166224 G > A, g.162166654 T > A, g.162167165 C > A, g.162167314 A > C, g.162172653 G > C and g.162182976 C > T). In block 2, the linkage disequilibrium of g.162182976 C > T with the other five SNPs (g.162166224 G > A, g.162166654 T > A, g.162167165 C > A, g.162167314 A > C, g.162172653 G > C) was not very strong (r 2 < 0.33).
In total, nine haplotypes were constructed. The haplotypes of the PRKG2 gene and their frequencies in the Dezhou donkey are shown in Table 3 were not used for association analysis. Hap2Hap6, Hap-2Hap8, Hap4Hap6, Hap4Hap7, Hap5Hap6, Hap5Hap7, Hap5Hap8, Hap5Hap9, Hap6Hap6, Hap7Hap8 and Hap-9Hap9 combinations were not found in our population.

Association analysis of PRKG2 SNPs with the number of thoracolumbar vertebrae and carcass traits in Dezhou donkeys
The association analysis of PRKG2 SNPs with the number of thoracolumbar vertebrae and carcass traits in Dezhou donkeys are shown in Table 4. The results of association analysis showed that the g.162149155 G > C and g.162158453 C > T mutations of the PRKG2 gene were significantly associated with the body height (P < 0.05). The g.162140112 A > G was significantly associated with differences in the number and length of lumbar vertebrae (P < 0.05). g.162153251 G > A (P < 0.01), g.162156524 C > T (P < 0.01), g.162158453 C > T (P < 0.05) and g.162163775 T > G (P < 0.05) were significantly associated with carcass weight. In addition to being significantly associated with body height and carcass weight, g.162158453 C > T was significantly associated with chest circumference (P < 0.05). Our analysis showed that there were significant relationships between the different locus of the g.162166224 G > A, g.162166654 T > A, g.162167165 C > A, g.162167314 A > C, g.162172653 G > C and the number of thoracic vertebrae in Dezhou donkey (P < 0.05). The g.162163775 T > G locus was significantly associated with the total number of thoracic and lumber, and the total number of thoracolumbar vertebrae was higher in donkeys with the TG genotype than in those with the TT genotype (P < 0.01).

Association analysis of PRKG2 haplotype combinations with the number of thoracolumbar vertebrae and carcass traits in Dezhou donkeys
Different haplotype combinations were not significantly associated with body height, body length, chest circumference, carcass weight, the number of lumbar vertebrae, the length of lumbar vertebrae, the number of thoracic vertebrae, the length of thoracic vertebrae, the total number of thoracic and lumbar vertebrae (P > 0.05) (

Discussion
The Dezhou donkey is one of China's five best donkey breeds, with high production characteristics and stable genetic performance [23]. In recent decades, breeding efforts have focused on animals that meet people's basic needs, such as pigs and chickens. After satisfying food and clothing, people's demand for food began to pursue nutrition and health. Many studies showed that donkey meat is of great nutritional value [24]. However, as a special-type economic animal, the progress of donkey breeding is slow. Therefore, the identification of molecular markers affecting economic traits is essential to accelerate the molecular breeding process of Dezhou donkeys. Fifteen SNPs were identified in the PRKG2 gene of the Dezhou donkey for the first time in this study, and SNPs located in the PRKG2 gene have not been previously reported in donkeys. Polymorphisms in the PRKG2 gene have also been found in humans, dogs, and cattle. The mutation c.1705 C > T found in the exonic region of the human genome is associated with acral dysplasia [25]. Koltesa et al. (2009) found that the C/T transition in exon 15 of the American Angus cattle PRKG2 gene introduced a stop codon (R678X) and demonstrated that the R678X resulted in the loss of regulation of COL2 and COL10 mRNA expression. R678X is a pathogenic mutation in American Angus cattle dwarfism. The fifteen SNPs we identified were all located in the intron region. Similarly, the c.1634 + 1 G > T locus found in the intron region of dogo Argentines is a candidate pathogenic variant of dwarfism. Radiographs of dogs with dwarfism show reduced levels of endochondral ossification in epiphyseal plates and premature closure of the distal ulna epiphysis line [11]. Currently, genetic variants of the PRKG2 gene have not been identified in horses, sheep, and pigs.
In the fifteen SNPs confirmed, only one mutant site was low polymorphic, six mutant sites were moderately polymorphic, and eight mutant sites were highly polymorphic. This result indicates a relatively high level of Values with different letters (a > b; A > B) within the same row denote significance levels of P < 0.05 and P < 0.01, respectively polymorphism in this population. However, considering that our group consisted entirely of two-year-old male donkeys, our results have some limitations. The g.162153251 G > A, g.162156524 C > T, and g.162216538 G > A locus are not in HWE, indicating that they may be affected by artificial selection, natural selection, migration, and population size, and the genetics of these three sites are unstable [26]. The average observed heterozygosity of fifteen SNPs was 0.3575, and the average expected heterozygosity was 0.6022, this suggests that the Dezhou donkey population is rich in genetic variation [27]. Growth traits are important indicators of breeding, and thirteen SNPs were significantly associated with the thoracolumbar number and carcass traits. Unfortunately, g.162160146 T > C and g.162182976 C > T were not associated with all traits; this may be due to the small sample size used in our study [19,21]. g.162149155 G > C and g.162158453 C > T were significantly associated with the body height of the Dezhou donkey (P < 0.05). Duyvenvoorde et al. (2014) showed that the PRKG2 gene was identified as a candidate gene for human height. However, fifteen SNPs of the PRKG2 gene were not significantly associated with body length in our study. g.162140112 A > G was significantly associated with lumbar spine number and length (P < 0.05). g.162166224 G > A, g.162166654 T > A, g.162167165 C > A, g.162167314 A > C and g.162172653 G > C were significantly associated with the number of thoracic vertebrae (P < 0.05), and g.162163775 T > G was significantly associated with the total number of thoracolumbar vertebrae (P < 0.01). Yi et al. (2021) found that the PRKG2 gene promotes adipogenesis and impairs osteoblastogenesis. It is the opposite of our results, g.162140112 A > G, g.162163775 T > G, g.162166224 G > A, g.162166654 T > A, g.162167165 C > A, g.162167314 A > C and g.162172653 G > C may affect the function of osteoclastogenesis in the PRKG2 gene has been hypothesized, but the mechanisms involved need to be further investigated.
Haplotype combinations are highly likely to be inherited together [26]. Although SNP sites were significantly associated with carcass traits and the number of thoracolumbar vertebrae, association analysis revealed that the constructed haplotype combinations were not significantly associated with the number of thoracolumbar vertebrae and carcass traits. A possible explanation for this is that haplotype combination with the highest value of traits had a sample size of less than 6 were not included in the association analysis of this study [28]. Furthermore, donkeys with haplotype combination Hap4Hap8 had the significant length of lumbar vertebrae, number of lumbar vertebrae, and the total number of thoracolumbar vertebrae compared to donkeys with other haplotype combinations. Donkeys with haplotype combination Hap3Hap5 had the greatest carcass weight, length of thoracic vertebrae, and the number of thoracic vertebrae compared to donkeys with other haplotype combinations. Although there were no significant differences between haplotype combinations and traits, the dominant haplotype combinations Hap4Hap8 and Hap3Hap5 that we found were able to bring about some positive effects.
Similarly, SNPs located in introns significantly associated with growth performance compared with SNPs located in exons and non-coding regions. For example, a novel g.3624 A > G polymorphism in intron 2 of the TBX3 gene is significantly associated with body size in donkeys [20]. Numerous studies have shown that SNPs located in introns are associated with alternative splicing. Alternative splicing plays a vital role in regulating biological functions [29]. The g.19970 A > G site found in intron 11 of the cow INCNEP gene enhances the action of the splicing factor SRSF1, SRSF1(IgM-BRCA1), and SRSF5. It changes the binding sites of splicing factor SRSF6, generating a new transcript that alters gene expression [30]. g.11043 C > T in the intron 1 of the SPEF2 gene that alters the binding of the splicing factor binding protein SC35 to the target sequence, and it was hypothesized that this mutation is essential for the production of new transcripts and therefore has an effect on bull semen trait production [31]. The fifteen SNPs that were newly identified by us affected the shear factor binding sites that need to be further confirmed.

Conclusions
In this study, we focused on the variation of the PRKG2 gene and its association with the number of thoracolumbar vertebrae and carcass traits of donkeys. Based on the targeted and Sanger sequencing methods, we found fifteen SNPs of the PRKG2 gene, all located in the intron region. The results showed that the PRKG2 gene could be a molecular marker with multiple thoracolumbar vertebrae and better carcass traits in donkeys, laying the foundation for breeding high-quality donkey breeds with high meat production.