INTRODUCTION
Growing demand for lean, nutritious meat has increased the popularity of goat meat worldwide [1]. In Thailand, the number of goats raised for meat production continues to increase, although most goats are owned by smallholder farmers [2]. The traditional goat production system in Thailand is extensive grazing. To meet the increasing demand, an intensive production system is needed; however, labor and financial constraints make it difficult for smallholder farmers to invest in equipment. Ahlawat et al [3] reported that increasing the reproductive rate of low-producing goat breeds could help bridge the gap between meat demand and production. Moreover, improved reproductive performance of livestock breeds increases the number of available offspring, which in turn increases farmers’ profitability. Thus, high fertility in goats is an important factor for farmers.
Litter size, defined as the total number of kids born per kidding per goat, is a key component trait of doe reproduction. It is positively correlated with parity in goats, and multiple birth data are used in selection [4]. Most traits related to reproduction, including litter size, have low heritability, and genetic improvement is slow when selection is based on one or a few phenotypic records [5]. Genetic marker-assisted selection based on molecular indicators may accelerate the genetic improvement of goat fecundity. Zhang et al [6] reported that 69 molecular markers have been identified to improve the litter size of goats; however, the presence and value of these markers vary among breeds. Moreover, de Lima et al [5] suggest that litter size traits involve multiple genes and loci incorporating environmental factors and maternal influences such as maternal age and the intrauterine environment. Therefore, genomic analysis with a large number of molecular markers across the entire genome, rather than single genes, is appropriate to find markers associated with doe reproduction traits.
The genotyping by sequencing (GBS) technique is a simple, highly multiplexed, reproducible, and cost-effective method for the simultaneous discovery and genotyping of many markers [7]. The application of the GBS method to livestock provides an adequate solution for population genetic and genomic studies. Regarding fertility traits, this method has been successfully used to identify selection signatures of litter size in sheep [8] and candidate genes in Indian buffalo [9]. In this study, we performed single nucleotide polymorphism (SNP) discovery using the GBS method and conducted an association analysis to identify SNPs and putative genes related to doe reproduction traits, including litter size and age at the first, second, and third parities, as well as the intervals between parities, in goats. Our results may help to further accelerate the genetic improvement of goats, and the obtained SNP markers can be used to select genetically superior individuals in goat production.
MATERIALS AND METHODS
Ethics statement
The study was approved by the Animal Care and Use Committee of the Suranaree University of Technology (SUT, Nakhon Ratchasima, Thailand) (SUT 4/2558). The research was conducted following regulations on animal experiments and the Guidelines for the Use of Animals in Research, as recommended by the National Research Council of Thailand (U1-08780-2563).
Animals
A total of 31 female goats with multiple breeding records (age and litter size of individuals in the first, second, and third parities) from the SUT Farm were used in this study. The 31 goats were born between 2012 and 2019 and reached their third parity at ages ranging from 2 years and 5 months to 6 years and 7 months. Of the 31 goats, there were one Anglo-Nubian purebred, one Boer purebred, and 29 crossbred goats. The crossbred goats were produced using a total of 7 breeds: 5 goats were crossed using 2 breeds, 22 goats were crossed using 3 breeds, and 2 goats were crossed using 4 breeds. The breeding information for each goat analyzed in this study is shown in Supplement 1. The breed ratios of most crossbred goats were different, with only two of 29 crossbred goats having the same breed ratio. The breed proportions of the 31 goats used in this study were: Anglo-Nubian 32%, Saanen 26%, Chami 18%, Boer 11%, Thai native breed 9%, Alpine 2%, and Kalahari Red 2%. All goats were born on the SUT farm and reared according to the farm’s animal husbandry standards for meat production. The farm uses natural or artificial insemination to breed only healthy females that are at least 8 months old and weigh at least 30 kg when in heat. The farm also uses breeding records up to the third parity for selection when culling female goats.
Sample collection and DNA extraction
Approximately 3.0 mL of blood was collected aseptically from the jugular vein of each doe into a blood collection tube containing K3 EDTA (Hebei Xinle Sci&Tech Co. Ltd., Hebei, China) in October 2022. Genomic DNA was extracted from blood using the DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany) according to the manufacturer’s instructions. The quality of DNA was measured by the NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Massachusetts, MA, USA). The integrity of DNA was checked by 1% agarose (w/v) gel electrophoresis.
Genotyping by sequencing
Novogene Biotechnology Company (Novogene, Beijing, China) conducted the DNA library preparation and subsequent sequencing. The 250 ng of genomic DNA from each sample was digested with 5 units of MseI in a 20 μL reaction. The reaction was incubated at 37°C for 2 hours, followed by heat inactivation of the enzyme at 65°C for 20 min. The resulting fragments were ligated with two barcode adapters that either with a compatible sticky end with the digestion enzyme and the Illumina P5 or P7 universal sequences. The ligation reaction was incubated at 22°C for 2 hours, followed by heat inactivation at 65°C for 20 min. After several rounds of polymerase chain reaction (PCR) amplification, all the samples were pooled, purified using magnetic beads, and size-selected. Library quality was assessed and quantified, and genome DNA libraries with an average insert size of ~350 bp were constructed. Libraries were sequenced on an Illumina Novaseq 6000 instrument (Illumina, San Diego, CA, USA) to generate 150 bp paired-end reads. The raw reads obtained from Illumina sequencing were quality-controlled to obtain clean reads. Specifically, adapter-containing read and paired reads with more than 10% uncertain nucleotides or more than 50% low-quality nucleotides (base quality less than 5) were discarded using fastp version 0.20.0 [10]. The clean reads were aligned to the reference genome of the goat, ARS1 (INSDC Assembly GCA_001704415.1), using Burrows-Wheeler Aligner version 0.7.8-r455 [11] with parameters: mem -t 4 -k 32 -M. SNPs were called using the command ‘mpileup’ in SAMtools version 0.1.19-44428cd [12] with the parameters set as ‘-m 2 -F 0.002 -d 1000’. SNPs with a depth of the variant position > 4 and mapping quality >20 were retained for future analysis.
Genotype data quality control
The genotype data quality was controlled using VCFtools [13] before the association analysis. An SNP was removed if it (i) had more than two alternative alleles, (ii) had a sequencing depth <10, (iii) had an assigned genotype quality <20, (iv) had more than 5% missing genotypes, or (v) had a minor allele frequency <0.05. After filtering by SNP characteristics, 21,665 SNPs were retained. The obtained dataset of 21,665 SNPs was imputed using the LD-kNNi method in Tassel version 5.2.94 [14] and used for subsequent analyses.
Population structure and linkage disequilibrium
With 21,665 SNPs, principal components (PCs), kinship, and linkage disequilibrium (LD) among the genome-wide markers were calculated using Tassel version 5.2.94 [14]. The same 21,665 SNPs were used in STRUCTURE [15], testing population subgroups with K = 1 to 5 to determine the optimal number of population subgroups. Using the admixture model in STRUCTURE, K = 1 to 5 were tested with 150,000 reps, declaring the first 50,000 as burn-ins and performing 10 iterations for each K value. To determine the optimal number of population subgroups, cross-validation analysis was performed using ADMIXTURE version 1.3.0 [16]. LD decay plot was generated based on the r2 values and the distance between each SNP pair.
Identification of single nucleotide polymorphisms associated with reproduction traits
All doe reproduction traits were tested for deviation from normality using the Shapiro–Wilk test. Before the association analysis, the square root transformation was used to analyze the non-normally distributed trait, the interval between the first and second parities in days. Litter size data, consisting of two levels (single and twin births), were treated as a qualitative trait. The litter size of individuals at the first, second, and third parities was converted with single birth coded as 0 and twin birth coded as 1. Association analysis was performed using the general linear model analysis in Tassel version 5.2.94 [14], and the first four PC vectors were used as covariates to account for population structure. The results were visualized in a Manhattan plot. The calculated p-values were subjected to false discovery rate (FDR) correction, and SNPs with FDR ≤ 0.05 were determined to be significant. For the significant SNPs, the goat reference genome (GCF_001704415.2) on the NCBI database was used to search for candidate genes within a 100 kb downstream and upstream window from the SNPs. The local LD was viewed for the significant SNPs, and LD plots were created using Haploview 4.2 [17]. The LD blocks were defined using the confidence interval method and default parameters. For litter size traits that showed correlation, genome-wide association analysis was performed using the R package rrBLUP [18], with parity as a fixed effect. The significance of SNPs was determined in the same manner as described above.
Genotyping method for candidate single nucleotide polymorphism
This work has targeted only SNP marker development. Genotypes of the significant SNP containing candidate genes were confirmed using the PCR allele competitive extension (PACE) genotyping assay. The PACE genotyping assay involves PCR with endpoint fluorescence measurement using two allele-specific forward primers and one common reverse primer [19]. Three primers were designed by 3CR Bioscience (Essex, UK): allele-specific primer 1 (HEX-GCACAGCTA GTATTTTGGTTTTCCC), allele-specific primer 2 (FAM-AGCACAGCTAGTATTTTGGTTTTCCA), and the common primer (GAGTGGGTTTGCCTGACGAGAAAA). The fluorescent PACE genotyping assay was performed on a CFX Opus Real-Time PCR System (Bio-Rad, Hercules, CA, USA) in a final volume of 10 μL, containing 1× PACE genotyping master mix (3CR Bioscience), 12 μM of a mix of extended allele-specific primers, 30 μM of common primer, and 1 to 10 ng of template DNA. The touch-down PCR amplification conditions for PACE genotyping included hot-start activation at 94°C for 15 min, followed by 10 cycles of 94°C for 20 s, 65°C for 60 s (dropping 0.8°C per cycle), then 30 cycles of 94°C for 20 s and 57°C for 60 s, with a final step at 37°C for 60 s. CFX Maestro Software 2.3 (Bio-Rad) was used to analyze SNP genotypes, distinguishing between two homozygous genotypes emitting either FAM or HEX and a heterozygous genotype emitting both fluorescences. The assay was conducted with 31 DNA samples and two negative controls.
RESULTS
Litter size and age for the three parities
The ages and numbers of kids born at the three parities of the 31 female goats are summarized in Table 1. Nine traits were considered as doe reproduction traits: litter size at the first, second, and third parities (LS1, LS2, and LS3), age in days at the first, second, and third parities (AP1, AP2, and AP3), and the intervals in days between the first and second, second and third, and first and third parities (I1, I2, and I3). The averages of LS1, LS2, and LS3 were 1.12, 1.45, and 1.77, with standard deviations of 0.34, 0.50, and 0.42, respectively. All goats gave birth to one or two kids during each of the three parities, and litter size tended to increase with parity. The averages of AP1, AP2, and AP3 were 762, 1156, and 1655 days, with standard deviations of 242, 313, and 356 days, respectively. The averages of I1, I2, and I3 were 393, 499, and 893 days, with standard deviations of 219, 276, and 309 days, respectively. The nine reproduction traits of the individual goats analyzed in this study are presented in Supplement 2. Correlation between doe reproduction traits was examined using Spearman’s rank correlation coefficient test (Figure 1). A moderately significant correlation was found between LS1 and LS2.
Sequencing, mapping, and single nucleotide polymorphism calling
Sequence quality assessment and SNP calling results are summarized in Table 2. After filtering out low-quality reads, an average of 9,263,800 reads were mapped to the goat reference genome, with a mapping rate of 99.72%. An average of 304,852 SNPs were initially detected in the 31 female goats. The sequencing, mapping, and SNP calling results for individual goats are presented in Supplement 3. The 21,665 SNPs were retained after filtering by SNP characteristics. The information about each SNP, including physical distance between SNPs, heterozygosity, and SNP distribution across each chromosome, is presented in Supplements 4, 5.
Genomic relatedness and population structure in studied goats
The 21,665 SNPs were used to examine the genetic structure of the 31 goats. The relatedness between 31 goat individuals was illustrated by PC analysis plots and kinship heatmap (Figures 2, 3). In the figure, each goat is classified and colored based on the breeds used in its production. The PC analysis revealed genetic differentiation, with the first two components explaining 19.90% of the genetic variation (11.39% for PC1 and 8.51% for PC2) (Figure 2). The purebred goats, Anglo-Nubian and Boer, were plotted at opposite positions on the PC1 axis (red and blue color, respectively). The crossbred goat that was bred with 87.5% Saanen and 12.5% Thai native breed was plotted in the lower left of the figure (dark green). Figure 3 facilitates understanding the genetic relationships between individuals. The 31 goats were divided into two large clusters based on kinship value. The first cluster contained 18 goats and the second cluster contained 13 goats, with kinship values of 0.02 and 0.08, respectively. The heatmap displays the self-relatedness along the diagonal and relatedness between individuals outside the diagonal, with each breed group as an axis. A high relatedness, kinship value higher than 0.25, was observed in 15 pairs, indicating the 31 female goats include individuals related by pedigree. The information about the PC of each goat and kinship is shown in Supplements 6, 7. The model parameters for determining the optimal number of populations obtained from STRUCTURE and ADMIXTURE are shown in Table 3. Results from STRUCTURE suggested K = 2 as the most likely scenario, however, the ADMIXTURE cross-validation analysis indicated that K = 1 was the optimal modeling choice, revealing a weak genetic structure of the 31 goats. The LD half-decay distance of 481 bp was found in the 31 goat population (Figure 4).
Detection of single nucleotide polymorphisms associated with doe reproduction traits
The results of the association analysis of the doe reproduction traits are shown in Figure 5. Among the nine traits, significant SNPs were detected only for the LS1 trait (Figure 5A). No SNPs were found to be significantly associated with LS2, LS3, AP1, AP2, AP3, I1, I2, and I3 (Figures 5B–5I). A total of six SNPs on chromosomes 4, 11, 24, and 29 exhibited a significant association with the LS1 trait (Table 4). Five SNPs identified genes within 100 kb downstream and upstream of the significant SNPs. The most significant SNP was detected on chromosome 4, with three genes identified around it: IKAROS family zinc finger 1 (IKZF1), fidgetin-like 1 (FIGNL1), and dopa decarboxylase (DDC). Three SNPs detected on chromosome 29 were located on the asparaginase and isoaspartyl peptidase 1 (ASRGL1) gene.
LD plots were generated for the significant SNPs on chromosomes 4, 11, 24, and 29 to examine the local LD blocks. No LD block containing significant SNP was identified on chromosomes 4, 11, and 24 (Figures 6, 7). Three significant SNPs on chromosome 29 were located within an LD block (Figure 8). However, no other SNPs were detected within the same LD block. Since a moderately significant correlation was observed between LS1 and LS2, further association analysis was performed by combining the LS1 and LS2 data as the phenotype and adjusting for parity as a fixed effect. Although no significant SNPs were identified, the most significant SNP (S4_115015047) in the LS1 trait analysis also had the lowest p-value in this analysis, slightly falling below the suggestive significance threshold (Figure 9).
Development of a genotyping method for candidate single nucleotide polymorphism
The genotypes at the most significant SNP (S4_115015047) for the LS1 trait were confirmed using the PACE genotyping assay (Figure 10). The PACE genotyping assay identified three genotypes: 26 individuals with genotype TT, 4 with genotype GT, and 1 with genotype GG. In contrast, the genotyping data obtained from the GBS method showed 25 individuals with genotype TT, 5 with genotype GT, and missing data for 1 individual. The genotypes determined by the GBS method and the PACE genotyping assay showed nearly identical results, with a 96% concordance rate (29/30) when excluding the missing genotype data (Supplement 8).
DISCUSSION
Improving goat reproductive performance enhances farmers’ productivity and profitability. Identifying DNA variants associated with goat reproductive traits will enable selection based on these variants and accelerate genetic improvement. In this study, the GBS library from 31 female goats successfully generated 304,852 SNPs. However, this method is limited by the high amount of missing data and insufficient coverage depth, which often results in incorrect genotyping. Minimizing missing data can be achieved by sequencing fewer sites with higher target coverage [20]. Consequently, we selected 21,665 SNPs by primarily filtering based on sequencing depth and then examined the association between these SNPs and nine doe reproduction traits.
The PC plot obtained using the 21,665 SNPs confirmed three clusters, explaining 19.90% of the variation when the first two components were combined. The purebred Boer goat is plotted where PC2 values are higher than 10, the purebred Anglo-Nubian goat is placed where PC2 values are between −10 and 10, and the crossbred goat with a high proportion of Saanen breed is located where PC2 values are lower than −10, suggesting that PC2 contains components related to breed characteristics. The two main clusters generated based on kinship values included purebred Anglo-Nubian and Boer breeds in their respective clusters. All goats in the cluster belonging to the purebred Anglo-Nubian breed had PC1 values greater than −6. The other goats in the cluster belonging to the purebred Boer breed had PC1 values less than −8, except for one individual, NG242_22. This suggests that the components of PC1 include kinship values. In a study of Tunisian sheep using 115,121 SNPs obtained through the GBS method, Bedhiaf-Romdhani et al [21] reported that their PCA results showed the first two components together explained 20.8% of the variation, revealing clear differentiation between two European-origin breeds and four other breeds. These findings suggest that even a relatively small number of SNPs can reveal genetic characteristics by breed in the female goats studied. Our population structure analysis determined K = 1, suggesting that the 31 goats form a single genetic group with no clear subpopulation differentiation. This indicated extensive admixture due to interbreeding of the seven breeds in various combinations. The rapid LD decay of 481 bp reflects high genetic diversity within the population, which is consistent with the unique breed ratios observed in most crossbred individuals. These results represent the genetic characteristics of the 31 female goats used in this study, specifically the genetic homogeneity at the population level and genetic diversity at the individual level.
Genome-wide association studies (GWAS) are a powerful tool for detecting the genetic basis of complex traits. Adequate statistical power for GWAS requires a larger sample size and appropriate case-control designs. In this study, 31 female goats with multiple birth records were used due to the importance of these traits for on-farm female goat selection. Previous studies have successfully identified candidate genes associated with litter size in goats using relatively small sample sizes (31 and 40) from single farms [22,23], suggesting that key genes may also be detectable in this study. Among the 21,665 SNPs analyzed, only six SNPs were associated with litter size at the first parity (LS1). Of these, three SNPs formed one LD block and were located on gene ASRGL1. ASRGL1 is an enzyme involved in L-aspartic acid production, and Fitzgerald et al [24] reported that ASRGL1 protein was specifically expressed in the proliferating uterine fluid of infertile women. However, little is known about its role in fertility.
Although no LD block was identified, the most significant SNP associated with LS1 was detected on chromosome 4 (S4_115015047), and three genes—IKZF1, FIGNL1, and DDC—were identified within 100 kbp of it. IKAROS (IKZF1) is a transcription factor and a member of the IKAROS gene family. Members of this family play an important role in lymphocyte development and proliferative responses [25] and are significant regulators in the hematopoietic system [26]. They regulate the expression of target genes by binding to promoter regions and play a role in the gene regulatory network in the human body [27]. In mouse pituitary tissue, IKAROS is expressed and regulates the expression of adrenocorticotropic hormone and the secretion of adrenal cortical hormones [28]. However, the role of IKAROS in other tissues is largely unknown. In goats, Zhang et al [29] suggested that IKZF1 may share similar functions with RUNX3 by analyzing single-nucleus transcriptomics data from the ovaries of polytocous and monotocous goats (Nubian and Du’an goats, respectively) during the follicular phase. In female mice, the RUNX3 transcription factor regulates ovarian function and ovulation [30]. Moreover, Ojima et al [30] demonstrated that RUNX3 regulates folliculogenesis and steroidogenesis in granulosa cells of immature mice. Follicle maturation and ovulation rate are reported to be critical factors affecting litter size, increasing fertilization rates and the likelihood of more offspring [31]. Fidgetin (FIGN) is an ATP-dependent microtubule-associated factor involved in chromosome segregation and cell division that promotes rapid reassembly and nucleation of microtubules from the centrosome [32]. FIGNL1 is a gene involved in homologous recombination repair, which in mice causes male infertility due to meiotic defects [33]. In female mice, analysis of in vitro mature oocytes by Li et al [34] showed that knockdown of FIGN results in multisperm fertilization. Moreover, they found that FIGNL1 complemented the role of FIGN in female FIGN knockout mice, resulting in healthy and normal fertility. In female rats, the expression of FIGNL1 protein in primordial follicles obtained from the ovarian tissue of aged rats was downregulated compared to those of immature rats [35]. In humans, variations in the FIGNL1 gene have been associated with premature ovarian insufficiency, a disorder characterized by the cessation of menstrual cycles before the age of 40 due to the depletion or dysfunction of ovarian follicles [36]. The DDC gene is responsible for the synthesis of dopamine, which is essential for a healthy pregnancy [37]. Gratz et al [37] found that patients with recurrent miscarriages had lower levels of DDC and the D2-dopamine receptor in the trophoblast and decidua compared to healthy controls. In an analysis comparing the transcriptomes of 12 brain regions in 1-month pregnant and non-pregnant Murciano-Granadina goats, Luigi-Sierra et al [38] found that DDC was downregulated in the brainstem and olfactory bulb of 1-month pregnant goats. Dopamine is also a key modulator of maternal behavior in rodents [39]. In Small Tail Han sheep, Wang et al [40] reported that the levels of neurotransmitters, including dopamine, in the serum of ewes with twin lambs were significantly higher than those of ewes with single lambs (p<0.05), and ewes that had twin lambs showed better maternal behavior quality. These findings support the notion that the three genes (IKZF1, FIGN1, and DDC) identified near the most significant SNP in this study may be associated with fertility in female goats. Unlike the other five SNPs, the SNP S4_115015047 also showed the lowest p-value in the association analysis considering the number of parities, suggesting this SNP plays a role in determining litter size between the first and second parities. However, the mechanism by which genotype differences at the detected SNP relate to goat litter size remains unclear.
Genotypes of the SNP S4_115015047 associated with LS1, as obtained by the GBS method, were confirmed for the same individuals by the PACE genotyping assay, yielding a 96% concordance rate. One discrepancy was observed where an individual was heterozygous for genotype GT by the GBS method but homozygous for genotype GG by the PACE genotyping assay. This discrepancy may be attributed to differences in methodology; the PACE genotyping is a targeted method using allele-specific primers and does not require the complex bioinformatics involved in GBS variant calling. Interestingly, all individuals determined to be heterozygous for genotype GT by the PACE genotyping assay had twins at their first parity in this study. To test the usefulness of this marker, it is necessary to check the genotypes of other individuals; however, this was not possible due to the availability of samples in this study. Increasing goat litter size through selection improves production efficiency and provides economic benefits to farmers. The identified SNP and the applied PACE genotyping methods in this study could potentially be used to select female goats for future genetic improvement programs. However, because the litter size data analyzed in this study consisted of two levels, association analysis applying logistic regression would be necessary for a rigorous analysis. Moreover, future studies are required to further elucidate the significance of the mutation discovery by examining the genotypes of other female goats for the SNPs significantly associated with the litter size in this study.
CONCLUSION
Litter size variability is a significant economic trait with low heritability in domestic animals. In this study, we performed an association analysis in 31 goats using 21,665 SNPs obtained through the GBS method, focusing on nine doe reproduction traits. We identified six SNPs significantly associated with litter size at first parity. Three candidate genes, IKZF1, FIGNL1, and DDC, were located within 100 kb of the most significant SNP, suggesting their potential biological relevance. Furthermore, a straightforward genotyping method with allele-specific primers was established for this SNP, and all individuals identified as heterozygous for genotype GT had twins at their first parity. These findings may advance marker-assisted selection in goat breeding programs. However, the study’s limitation includes the relatively small sample size, attributed to the challenges in obtaining female goats with extensive breeding records and the limited functional understanding of the detected genes in goats. Therefore, further study on the frequency of genetic mutation with larger sample sizes and functional studies of candidate genes are necessary to further elucidate the significance of the mutation discovery and understand the mechanisms by which genotype differences affect goat fertility.
PDF Links
PubReader
ePub Link
Full text via DOI
Full text via PMC
Download Citation
Supplement1
Print
